From da707a78edfbb8d1089d7d98a53ffbd33a85f8ab Mon Sep 17 00:00:00 2001 From: Nicolas Poilvert Date: Sun, 12 Feb 2012 00:09:49 -0500 Subject: [PATCH 01/20] 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' From b4a1bf18ba5c09ee58801e2c373c3db1710b09b0 Mon Sep 17 00:00:00 2001 From: Nicolas Poilvert Date: Sun, 12 Feb 2012 00:10:40 -0500 Subject: [PATCH 02/20] remove empty line in __init__ --- skimage/util/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/skimage/util/__init__.py b/skimage/util/__init__.py index aa8ead55..d01470bd 100644 --- a/skimage/util/__init__.py +++ b/skimage/util/__init__.py @@ -1,3 +1,2 @@ from .dtype import * from .array_views import * - From 39a487834c0dbb4d9ee034573a887d3e93e94868 Mon Sep 17 00:00:00 2001 From: Nicolas Poilvert Date: Sun, 12 Feb 2012 01:06:47 -0500 Subject: [PATCH 03/20] more tests on rolling views --- skimage/util/array_views.py | 9 ++- skimage/util/tests/test_array_views.py | 76 +++++++++++++++++++++++++- 2 files changed, 81 insertions(+), 4 deletions(-) diff --git a/skimage/util/array_views.py b/skimage/util/array_views.py index 26c5a191..95467c94 100644 --- a/skimage/util/array_views.py +++ b/skimage/util/array_views.py @@ -209,9 +209,12 @@ def rolling_view(arr, window_shape): """ # -- basic requirements on inputs - assert isinstance(arr, np.ndarray) - assert isinstance(window_shape, tuple) - assert len(window_shape) == arr.ndim + if not isinstance(arr, np.ndarray): + raise ValueError('the input should be an ndarray object') + if not isinstance(window_shape, tuple): + raise ValueError('the window shape should be a tuple') + if not (len(window_shape) == arr.ndim): + raise ValueError('array dimension and window length dont match') # -- input array dimension N = arr.ndim diff --git a/skimage/util/tests/test_array_views.py b/skimage/util/tests/test_array_views.py index 0d062f59..913224c4 100644 --- a/skimage/util/tests/test_array_views.py +++ b/skimage/util/tests/test_array_views.py @@ -1,7 +1,7 @@ import numpy as np from nose.tools import raises from numpy.testing import assert_equal -from skimage.util.array_views import block_view +from skimage.util.array_views import block_view, rolling_view @raises(ValueError) @@ -65,3 +65,77 @@ def test_block_view_3D_array(): [58, 59]]], [[[76, 77], [82, 83]]]])) + + +@raises(ValueError) +def test_rolling_view_input_not_array(): + + A = [1, 2, 3, 4, 5] + rolling_view(A, (2,)) + + +@raises(ValueError) +def test_rolling_view_window_not_tuple(): + + A = np.arange(10) + rolling_view(A, [2]) + + +@raises(ValueError) +def test_rolling_view_wrong_window_dimension(): + + A = np.arange(10) + rolling_view(A, (2,2)) + + +@raises(ValueError) +def test_rolling_view_negative_window_length(): + + A = np.arange(10) + rolling_view(A, (-1,)) + + +@raises(ValueError) +def test_rolling_view_window_too_large(): + + A = np.arange(10) + rolling_view(A, (11,)) + + +def test_rolling_view_1D(): + + A = np.arange(10) + window_shape = (3,) + B = rolling_view(A, window_shape) + assert_equal(B, np.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]])) + + +def test_rolling_view_2D(): + + A = np.arange(5*4).reshape(5, 4) + window_shape = (4, 3) + B = rolling_view(A, window_shape) + assert_equal(B.shape, (2, 2, 4, 3)) + assert_equal(B, np.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]]]])) From 41a259ee26ee3d7e174d41fcb91dd5a148548b70 Mon Sep 17 00:00:00 2001 From: Nicolas Poilvert Date: Sun, 12 Feb 2012 12:24:09 -0500 Subject: [PATCH 04/20] simplifying code for clarity and modifying tests --- skimage/util/array_views.py | 99 +++++++++++--------------- skimage/util/tests/test_array_views.py | 8 +-- 2 files changed, 47 insertions(+), 60 deletions(-) diff --git a/skimage/util/array_views.py b/skimage/util/array_views.py index 95467c94..35a1286c 100644 --- a/skimage/util/array_views.py +++ b/skimage/util/array_views.py @@ -9,21 +9,17 @@ 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. + """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 + 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 + each element in the tuple represents the number of input array elements + to include in a block along the corresponding direction Returns ------- @@ -32,6 +28,7 @@ def block_view(arr, block): Examples -------- >>> import numpy as np + >>> import block_view >>> A = np.arange(4*4).reshape(4,4) >>> A array([[ 0, 1, 2, 3], @@ -76,14 +73,10 @@ def block_view(arr, block): [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') + raise TypeError('block needs to be a tuple') # -- basic invalid values for 'block' block_shape = np.array(block).astype(np.int) @@ -114,21 +107,18 @@ def block_view(arr, block): 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. +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_shape: N-tuple - tuple of size N that gives the shape of the elementary - window + window: N-tuple + tuple of size N that gives the shape of the elementary window Returns ------- @@ -136,21 +126,19 @@ def rolling_view(arr, window_shape): 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. + 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. + 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 ---------- @@ -158,11 +146,13 @@ def rolling_view(arr, window_shape): Examples -------- + >>> import numpy as np + >>> import rolling_view >>> 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) + >>> window = (3,) + >>> B = rolling_view(A, window) >>> B.shape (8, 3) >>> B @@ -181,8 +171,8 @@ def rolling_view(arr, window_shape): [ 8, 9, 10, 11], [12, 13, 14, 15], [16, 17, 18, 19]]) - >>> window_shape = (4, 3) - >>> B = rolling_view(A, window_shape) + >>> window = (4, 3) + >>> B = rolling_view(A, window) >>> B.shape (2, 2, 4, 3) >>> B @@ -210,28 +200,25 @@ def rolling_view(arr, window_shape): # -- basic requirements on inputs if not isinstance(arr, np.ndarray): - raise ValueError('the input should be an ndarray object') - if not isinstance(window_shape, tuple): - raise ValueError('the window shape should be a tuple') - if not (len(window_shape) == arr.ndim): + 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') - # -- input array dimension - N = arr.ndim + # -- defining some variables + arr_shape = np.array(arr.shape) + window_shape = np.array(window, dtype=arr_shape.dtype) # -- 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 ((arr_shape - window_shape) < 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') + if ((window_shape - 1) < 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 + out_shape = tuple(arr_shape - window_shape + 1) + window # -- strides of output 'rolling view' array out_strides = arr.strides + arr.strides diff --git a/skimage/util/tests/test_array_views.py b/skimage/util/tests/test_array_views.py index 913224c4..c0f97b80 100644 --- a/skimage/util/tests/test_array_views.py +++ b/skimage/util/tests/test_array_views.py @@ -4,7 +4,7 @@ from numpy.testing import assert_equal from skimage.util.array_views import block_view, rolling_view -@raises(ValueError) +@raises(TypeError) def test_block_view_block_not_a_tuple(): A = np.arange(10) @@ -15,7 +15,7 @@ def test_block_view_block_not_a_tuple(): def test_block_view_negative_shape(): A = np.arange(10) - block_view(A, (-2)) + block_view(A, (-2,)) @raises(ValueError) @@ -67,14 +67,14 @@ def test_block_view_3D_array(): [82, 83]]]])) -@raises(ValueError) +@raises(TypeError) def test_rolling_view_input_not_array(): A = [1, 2, 3, 4, 5] rolling_view(A, (2,)) -@raises(ValueError) +@raises(TypeError) def test_rolling_view_window_not_tuple(): A = np.arange(10) From 93cd71553e44a74be5f368c3ecc1b10a29db2d3e Mon Sep 17 00:00:00 2001 From: Nicolas Poilvert Date: Sun, 12 Feb 2012 12:52:01 -0500 Subject: [PATCH 05/20] make doctests pass with nose --- skimage/util/array_views.py | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/skimage/util/array_views.py b/skimage/util/array_views.py index 35a1286c..be9072da 100644 --- a/skimage/util/array_views.py +++ b/skimage/util/array_views.py @@ -28,7 +28,7 @@ def block_view(arr, block): Examples -------- >>> import numpy as np - >>> import block_view + >>> from skimage.util.array_views import block_view >>> A = np.arange(4*4).reshape(4,4) >>> A array([[ 0, 1, 2, 3], @@ -47,28 +47,29 @@ def block_view(arr, block): [ 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) + (4, 2, 3, 1, 2, 2) >>> B[2:, 0, 2] array([[[[52, 53], [58, 59]]], - + + [[[76, 77], [82, 83]]]]) """ @@ -87,6 +88,8 @@ def block_view(arr, block): 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) @@ -147,7 +150,7 @@ def rolling_view(arr, window): Examples -------- >>> import numpy as np - >>> import rolling_view + >>> from skimage.util.array_views import rolling_view >>> A = np.arange(10) >>> A array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) @@ -180,18 +183,18 @@ def rolling_view(arr, window): [ 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], @@ -206,6 +209,8 @@ def rolling_view(arr, window): 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) From 1108bd892e83a1fa6ea340e98587abaddd7cc10f Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Sun, 12 Feb 2012 14:59:06 -0500 Subject: [PATCH 06/20] 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], From c1c954e19d4d6e56ecdd405a1b58ad073688b6dc Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Sun, 12 Feb 2012 14:59:58 -0500 Subject: [PATCH 07/20] STY: pep8 in util.shape + its tests --- skimage/util/tests/test_shape.py | 36 ++++++++++++++++---------------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/skimage/util/tests/test_shape.py b/skimage/util/tests/test_shape.py index 127487f1..3b5f8d41 100644 --- a/skimage/util/tests/test_shape.py +++ b/skimage/util/tests/test_shape.py @@ -29,7 +29,7 @@ def test_view_as_blocks_block_too_large(): def test_view_as_blocks_wrong_block_dimension(): A = np.arange(10) - view_as_blocks(A, (2,2)) + view_as_blocks(A, (2, 2)) @raises(ValueError) @@ -49,17 +49,17 @@ def test_view_as_blocks_1D_array(): def test_view_as_blocks_2D_array(): - A = np.arange(4*4).reshape(4,4) - B = view_as_blocks(A, (2,2)) - assert_equal(B[0,1], np.array([[2, 3], + A = np.arange(4 * 4).reshape(4, 4) + 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_view_as_blocks_3D_array(): - A = np.arange(4*4*6).reshape(4,4,6) - B = view_as_blocks(A, (1,2,2)) + A = np.arange(4 * 4 * 6).reshape(4, 4, 6) + 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]]], @@ -85,7 +85,7 @@ def test_view_as_windows_window_not_tuple(): def test_view_as_windows_wrong_window_dimension(): A = np.arange(10) - view_as_windows(A, (2,2)) + view_as_windows(A, (2, 2)) @raises(ValueError) @@ -119,23 +119,23 @@ def test_view_as_windows_1D(): def test_view_as_windows_2D(): - A = np.arange(5*4).reshape(5, 4) + A = np.arange(5 * 4).reshape(5, 4) window_shape = (4, 3) 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], - [ 8, 9, 10], + assert_equal(B, np.array([[[[0, 1, 2], + [4, 5, 6], + [8, 9, 10], [12, 13, 14]], - [[ 1, 2, 3], - [ 5, 6, 7], - [ 9, 10, 11], + [[1, 2, 3], + [5, 6, 7], + [9, 10, 11], [13, 14, 15]]], - [[[ 4, 5, 6], - [ 8, 9, 10], + [[[4, 5, 6], + [8, 9, 10], [12, 13, 14], [16, 17, 18]], - [[ 5, 6, 7], - [ 9, 10, 11], + [[5, 6, 7], + [9, 10, 11], [13, 14, 15], [17, 18, 19]]]])) From 96f5d14000526727b453b50abfd5361a662a9ee7 Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Sun, 12 Feb 2012 15:06:58 -0500 Subject: [PATCH 08/20] FIX: revert version.py --- skimage/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/version.py b/skimage/version.py index d9d6b43c..ab15b0be 100644 --- a/skimage/version.py +++ b/skimage/version.py @@ -1,2 +1,2 @@ # THIS FILE IS GENERATED FROM THE SKIMAGE SETUP.PY -version='0.5dev' +version='unbuilt-dev' From cae47eda18efec41a077b2b37ff515dbf4db8898 Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Mon, 13 Feb 2012 09:37:44 -0500 Subject: [PATCH 09/20] STY: fix typos --- skimage/util/shape.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/util/shape.py b/skimage/util/shape.py index 22156e06..c0eeb40d 100644 --- a/skimage/util/shape.py +++ b/skimage/util/shape.py @@ -9,7 +9,7 @@ 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). + """Block view of the input n-dimensional array (using re-striding). Parameters ---------- @@ -112,7 +112,7 @@ def view_as_windows(arr_in, window_shape): window_shape: tuple Defines the shape of the elementary n-dimensional orthotope - (better know as hyperrectangle [1]) of the rolling window view. + (better know as hyperrectangle [1]_) of the rolling window view. Returns ------- From 74864966ac10fe730fd38372a57dd99ce7461ec1 Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Mon, 13 Feb 2012 09:41:02 -0500 Subject: [PATCH 10/20] STY: move 'util.shape' contributors to CONTRIBUTORS.txt --- CONTRIBUTORS.txt | 4 ++++ skimage/util/shape.py | 4 ---- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 8f479945..bea4d230 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -3,6 +3,7 @@ - Nicolas Pinto Colour spaces and filters + Shape views: ``shape.view_as_windows`` and ``shape.view_as_blocks`` - Damian Eads Morphological operators @@ -93,3 +94,6 @@ - Gaël Varoquaux Harris corner detector + +- Nicolas Poilvert + Shape views: ``shape.view_as_windows`` and ``shape.view_as_blocks`` diff --git a/skimage/util/shape.py b/skimage/util/shape.py index c0eeb40d..7d796405 100644 --- a/skimage/util/shape.py +++ b/skimage/util/shape.py @@ -1,7 +1,3 @@ -# Authors: Nicolas Poilvert -# Nicolas Pinto -# License: BSD 3-clause - __all__ = ['view_as_blocks', 'view_as_windows'] import numpy as np From 99efa1b771bb5679737e93093dee9aaea42b6899 Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Mon, 13 Feb 2012 21:20:31 -0500 Subject: [PATCH 11/20] ENH: add montage2d with tests --- skimage/util/montage.py | 87 ++++++++++++++++++++++++++++++ skimage/util/tests/test_montage.py | 59 ++++++++++++++++++++ 2 files changed, 146 insertions(+) create mode 100644 skimage/util/montage.py create mode 100644 skimage/util/tests/test_montage.py diff --git a/skimage/util/montage.py b/skimage/util/montage.py new file mode 100644 index 00000000..15646440 --- /dev/null +++ b/skimage/util/montage.py @@ -0,0 +1,87 @@ +__all__ = ['montage2d'] + +import numpy as np + + +def montage2d(arr_in, fill='mean'): + """Create a 2-dimensional 'montage' from a 3-dimensional input array + representing an ensemble of equally shaped 2-dimensional images. + + For example, montage2d(arr_in, fill) with the following `arr_in` + + +---+---+---+ + | 1 | 2 | 3 | + +---+---+---+ + + will return: + + +---+---+ + | 1 | 2 | + +---+---+ + | 3 | * | + +---+---+ + + Where the '*' patch will be determined by the `fill` parameter. + + Parameters + ---------- + arr_in: ndarray, shape=[n_images, height, width] + 3-dimensional input array representing an ensemble of n_images + of equal shape (i.e. [height, width]). + + fill: float or 'mean' + How to fill the 2-dimensional output array when sqrt(n_images) + is not an integer. If 'mean' is chosen, then fill = arr_in.mean(). + + Returns + ------- + arr_out: ndarray, shape=[alpha * height, alpha * width] + Output array where 'alpha' has been determined automatically to + fit (at least) the `n_images` in `arr_in`. + + Example + ------- + >>> import numpy as np + >>> from skimage.util.montage import montage2d + >>> arr_in = np.arange(3 * 2 * 2).reshape(3, 2, 2) + >>> print arr_in + [[[ 0 1] + [ 2 3]] + + [[ 4 5] + [ 6 7]] + + [[ 8 9] + [10 11]]] + >>> arr_out = montage2d(arr_in) + >>> print arr_out.shape + (4, 4) + >>> print arr_out + [[ 0. 1. 4. 5. ] + [ 2. 3. 6. 7. ] + [ 8. 9. 5.5 5.5] + [ 10. 11. 5.5 5.5]] + >>> print arr_in.mean() + 5.5 + """ + assert arr_in.ndim == 3 + + n_images, height, width = arr_in.shape + + # -- determine alpha + alpha = int(np.ceil(np.sqrt(n_images))) + + # -- fill missing patches + if fill == 'mean': + fill = arr_in.mean() + + n_missing = int((alpha ** 2.) - n_images) + missing = np.ones((n_missing, height, width), dtype=arr_in.dtype) * fill + arr_out = np.vstack((arr_in, missing)) + + # -- reshape to 2d montage, step by step + arr_out = arr_out.reshape(alpha, alpha, height, width) + arr_out = arr_out.swapaxes(1, 2) + arr_out = arr_out.reshape(alpha * height, alpha * width) + + return arr_out diff --git a/skimage/util/tests/test_montage.py b/skimage/util/tests/test_montage.py new file mode 100644 index 00000000..cd7151ab --- /dev/null +++ b/skimage/util/tests/test_montage.py @@ -0,0 +1,59 @@ +from nose.tools import assert_equal, raises +from numpy.testing import assert_array_equal + +import numpy as np +from skimage.util.montage import montage2d + + +def test_simple(): + n_images = 3 + height, width = 2, 3, + arr_in = np.arange(n_images * height * width) + arr_in = arr_in.reshape(n_images, height, width) + + arr_out = montage2d(arr_in) + + gt = np.array( + [[ 0. , 1. , 2. , 6. , 7. , 8. ], + [ 3. , 4. , 5. , 9. , 10. , 11. ], + [ 12. , 13. , 14. , 8.5, 8.5, 8.5], + [ 15. , 16. , 17. , 8.5, 8.5, 8.5]] + ) + + assert_array_equal(arr_out, gt) + + +def test_fill(): + n_images = 3 + height, width = 2, 3, + arr_in = np.arange(n_images * height * width) + arr_in = arr_in.reshape(n_images, height, width) + + arr_out = montage2d(arr_in, fill=0) + + gt = np.array( + [[ 0. , 1. , 2. , 6. , 7. , 8. ], + [ 3. , 4. , 5. , 9. , 10. , 11. ], + [ 12. , 13. , 14. , 0. , 0. , 0. ], + [ 15. , 16. , 17. , 0. , 0. , 0. ]] + ) + + assert_array_equal(arr_out, gt) + + +def test_shape(): + n_images = 15 + height, width = 11, 7 + arr_in = np.arange(n_images * height * width) + arr_in = arr_in.reshape(n_images, height, width) + + alpha = int(np.ceil(np.sqrt(n_images))) + + arr_out = montage2d(arr_in) + assert_equal(arr_out.shape, (alpha * height, alpha * width)) + + +@raises(AssertionError) +def test_error_ndim(): + arr_error = np.random.randn(1, 2, 3, 4) + montage2d(arr_error) From 64177b2000ea7a41a5a35d263ce00870f3693d79 Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Tue, 14 Feb 2012 00:04:42 -0500 Subject: [PATCH 12/20] ENH: add a 'normalize' kwarg to util.montage2d() --- CONTRIBUTORS.txt | 7 ++++--- skimage/util/montage.py | 18 ++++++++++++++++-- 2 files changed, 20 insertions(+), 5 deletions(-) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index bea4d230..5082c508 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -2,8 +2,9 @@ Project coordination - Nicolas Pinto - Colour spaces and filters - Shape views: ``shape.view_as_windows`` and ``shape.view_as_blocks`` + Colour spaces and filters. + Shape views: ``util.shape.view_as_windows`` and ``util.shape.view_as_blocks`` + Montage helpers: ``util.montage`` - Damian Eads Morphological operators @@ -96,4 +97,4 @@ Harris corner detector - Nicolas Poilvert - Shape views: ``shape.view_as_windows`` and ``shape.view_as_blocks`` + Shape views: ``util.shape.view_as_windows`` and ``util.shape.view_as_blocks`` diff --git a/skimage/util/montage.py b/skimage/util/montage.py index 15646440..0319fe84 100644 --- a/skimage/util/montage.py +++ b/skimage/util/montage.py @@ -2,8 +2,10 @@ __all__ = ['montage2d'] import numpy as np +EPSILON = 1e-6 -def montage2d(arr_in, fill='mean'): + +def montage2d(arr_in, fill='mean', normalize=True): """Create a 2-dimensional 'montage' from a 3-dimensional input array representing an ensemble of equally shaped 2-dimensional images. @@ -29,10 +31,13 @@ def montage2d(arr_in, fill='mean'): 3-dimensional input array representing an ensemble of n_images of equal shape (i.e. [height, width]). - fill: float or 'mean' + fill: float or 'mean', optional How to fill the 2-dimensional output array when sqrt(n_images) is not an integer. If 'mean' is chosen, then fill = arr_in.mean(). + normalize: bool, optional + Whether to normalize each image with zero-mean, unit-variance. + Returns ------- arr_out: ndarray, shape=[alpha * height, alpha * width] @@ -68,6 +73,15 @@ def montage2d(arr_in, fill='mean'): n_images, height, width = arr_in.shape + # -- normalize if necessary + if normalize: + arr_in = arr_in.T + arr_in -= arr_in.mean(0) + astd = arr_in.std(0) + astd[astd < EPSILON] = 1 + arr_in /= astd + arr_in = arr_in.T + # -- determine alpha alpha = int(np.ceil(np.sqrt(n_images))) From 96c2b3fc9d2039e0f4aa2ae5b7c23a0a0aad9650 Mon Sep 17 00:00:00 2001 From: Nicolas Poilvert Date: Mon, 13 Feb 2012 22:30:17 -0500 Subject: [PATCH 13/20] DOC: add example for util.shape.view_as_{blocks,windows} --- doc/examples/plot_view_as_blocks.py | 63 +++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) create mode 100644 doc/examples/plot_view_as_blocks.py diff --git a/doc/examples/plot_view_as_blocks.py b/doc/examples/plot_view_as_blocks.py new file mode 100644 index 00000000..305192dc --- /dev/null +++ b/doc/examples/plot_view_as_blocks.py @@ -0,0 +1,63 @@ +""" +============================ +Block views on images/arrays +============================ + +This example illustrates the use of `view_as_blocks` from `skimage.util.shape`. +Block views can be incredibly useful when one wants to perform local operations +on non-overlapping image patches. + +We use `lena` from `scipy.misc` and virtually 'slice' it into square blocks. +Then, on each block, we either pool the mean, the max or the median value of +that block. The results are displayed altogether, along with a 'classic' +`bicubic` rescaling of the original `lena` image. +""" + +import numpy as np +from scipy.misc import lena, imresize +from matplotlib import pyplot as plt +import matplotlib.cm as cm +from skimage.util.shape import view_as_blocks + + +# -- get `lena` from scipy in grayscale +l = lena() + +# -- size of blocks +block_shape = (4, 4) + +# -- see `lena` as a matrix of blocks (of shape +# `block_shape`) +view = view_as_blocks(l, block_shape) + +# -- collapse the last two dimensions in one +flatten_view = view.reshape(view.shape[0], view.shape[1], -1) + +# -- resampling `lena` by taking either the `mean`, +# the `max` or the `median` value of each blocks. +mean_view = np.mean(flatten_view, axis=2) +max_view = np.max(flatten_view, axis=2) +median_view = np.median(flatten_view, axis=2) + +# -- display resampled images +plt.figure(figsize=(10, 10)) + +plt.subplot(221) +plt.title("Original rescaled\n in bicubic mode"); +l_resized = imresize(l, view.shape[:2], interp='bicubic') +plt.imshow(l_resized, cmap=cm.Greys_r) + +plt.subplot(222) +plt.title("Block view with\n local mean pooling") +plt.imshow(mean_view, cmap=cm.Greys_r) + +plt.subplot(223) +plt.title("Block view with\n local max pooling"); +plt.imshow(max_view, cmap=cm.Greys_r) + +plt.subplot(224) +plt.title("Block view with\n local median pooling") +plt.imshow(median_view, cmap=cm.Greys_r) + +plt.subplots_adjust(hspace=0.4, wspace=0.4) +plt.show() From c715adf1c0cd8363df90d7f445cc105e795d7d62 Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Tue, 14 Feb 2012 00:16:00 -0500 Subject: [PATCH 14/20] DOC: add example of use of view_as_blocks with kmeans for simple dictionary learning --- doc/examples/plot_gabors_from_lena.py | 100 ++++++++++++++++++++++++++ 1 file changed, 100 insertions(+) create mode 100644 doc/examples/plot_gabors_from_lena.py diff --git a/doc/examples/plot_gabors_from_lena.py b/doc/examples/plot_gabors_from_lena.py new file mode 100644 index 00000000..0f483207 --- /dev/null +++ b/doc/examples/plot_gabors_from_lena.py @@ -0,0 +1,100 @@ +""" +======================================================= +Gabors / Primary Visual Cortex "Simple Cells" from Lena +======================================================= + +(under construction) + +How to build a (bio-plausible) "sparse" dictionary (or 'codebook', or +'filterbank') for e.g. image classification without any fancy math and +with just standard python scientific librairies? + +Please find below a short answer ;-) + +This simple example shows how to get Gabor-like filters [1]_ using just +the famous Lena image. Gabor filters are good approximations of the +"Simple Cells" [2]_ receptive fields [3]_ found in the mammalian primary +visual cortex (V1) (for details, see e.g. the Nobel-prize winning work of Hubel +& Wiesel done in the 60s). + +Here we use McQueen's 'kmeans' algorithm [4]_, as a simple bio-plausible +hebbian-like learning rule and we apply it (a) to patches of the +original Lena image (retinal projection), and (b) to patches of an +LGN-like [5]_ Lena image using a simple difference of gaussians (DoG) +approximation. + +Enjoy ;-) And keep in mind that getting Gabors on natural image patches +is not rocket science. + +.. [1] http://en.wikipedia.org/wiki/Gabor_filter +.. [2] http://en.wikipedia.org/wiki/Simple_cell +.. [3] http://en.wikipedia.org/wiki/Receptive_field +.. [4] http://en.wikipedia.org/wiki/K-means_clustering +.. [5] http://en.wikipedia.org/wiki/Lateral_geniculate_nucleus + +References +---------- +D. H. Hubel and T. N. Wiesel Receptive Fields of Single Neurones in the +Cat's Striate Cortex J. Physiol. pp. 574-591 (148) 1959 + +D. H. Hubel and T. N. Wiesel Receptive Fields, Binocular Interaction and +Functional Architecture in the Cat's Visual Cortex J. Physiol. 160 pp. +106-154 1962 +""" + +import numpy as np +from scipy import misc +from scipy.cluster.vq import kmeans2 +import matplotlib.pyplot as plt + +from skimage.util.shape import view_as_windows +from skimage.util.montage import montage2d +from scipy import ndimage as ndi + +np.random.seed(42) + +patch_shape = 8, 8 +n_filters = 49 + +lena = misc.lena() / 255. + +# -- filterbank1 on original Lena +patches1 = view_as_windows(lena, patch_shape) +patches1 = patches1.reshape(-1, patch_shape[0] * patch_shape[1])[::8] +fb1, _ = kmeans2(patches1, n_filters, minit='points') +fb1 = fb1.reshape((-1,) + patch_shape) +fb1_montage = montage2d(fb1) + +# -- filterbank2 LGN-like Lena +lena_dog = ndi.gaussian_filter(lena, .5) - ndi.gaussian_filter(lena, 1) +patches2 = view_as_windows(lena_dog, patch_shape) +patches2 = patches2.reshape(-1, patch_shape[0] * patch_shape[1])[::8] +fb2, _ = kmeans2(patches2, n_filters, minit='points') +fb2 = fb2.reshape((-1,) + patch_shape) +fb2_montage = montage2d(fb2) + +# -- +plt.figure(figsize=(9, 3)) + + +plt.subplot(2, 2, 1) +plt.imshow(lena, cmap=plt.cm.gray) +plt.axis('off') +plt.title("Lena (original)") + +plt.subplot(2, 2, 2) +plt.imshow(fb1_montage, cmap=plt.cm.gray) +plt.axis('off') +plt.title("K-means filterbank (codebook) on Lena (original)") + +plt.subplot(2, 2, 3) +plt.imshow(lena_dog, cmap=plt.cm.gray) +plt.axis('off') +plt.title("Lena (LGN-like DoG)") + +plt.subplot(2, 2, 4) +plt.imshow(fb2_montage, cmap=plt.cm.gray) +plt.axis('off') +plt.title("K-means filterbank (codebook) on Lena (LGN-like DoG)") + +plt.show() From 68f8c68c3ff286f60e07d571068e01fb70377483 Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Tue, 14 Feb 2012 19:36:07 -0500 Subject: [PATCH 15/20] FIX: util.montage2d 'normalize' becomes 'rescale_intensity' and defaults to False --- skimage/util/montage.py | 20 +++++++++----------- skimage/util/tests/test_montage.py | 22 ++++++++++++++++++++++ 2 files changed, 31 insertions(+), 11 deletions(-) diff --git a/skimage/util/montage.py b/skimage/util/montage.py index 0319fe84..a879a0a2 100644 --- a/skimage/util/montage.py +++ b/skimage/util/montage.py @@ -1,11 +1,12 @@ __all__ = ['montage2d'] import numpy as np +from .. import exposure EPSILON = 1e-6 -def montage2d(arr_in, fill='mean', normalize=True): +def montage2d(arr_in, fill='mean', rescale_intensity=False): """Create a 2-dimensional 'montage' from a 3-dimensional input array representing an ensemble of equally shaped 2-dimensional images. @@ -35,8 +36,8 @@ def montage2d(arr_in, fill='mean', normalize=True): How to fill the 2-dimensional output array when sqrt(n_images) is not an integer. If 'mean' is chosen, then fill = arr_in.mean(). - normalize: bool, optional - Whether to normalize each image with zero-mean, unit-variance. + rescale_intensity: bool, optional + Whether to rescale the intensity of each image to [0, 1]. Returns ------- @@ -73,14 +74,11 @@ def montage2d(arr_in, fill='mean', normalize=True): n_images, height, width = arr_in.shape - # -- normalize if necessary - if normalize: - arr_in = arr_in.T - arr_in -= arr_in.mean(0) - astd = arr_in.std(0) - astd[astd < EPSILON] = 1 - arr_in /= astd - arr_in = arr_in.T + # -- rescale intensity if necessary + if rescale_intensity: + for i in xrange(n_images): + arr_in[i] = exposure.rescale_intensity( + arr_in[i], out_range=(0.0, 1.0)) # -- determine alpha alpha = int(np.ceil(np.sqrt(n_images))) diff --git a/skimage/util/tests/test_montage.py b/skimage/util/tests/test_montage.py index cd7151ab..47e5426d 100644 --- a/skimage/util/tests/test_montage.py +++ b/skimage/util/tests/test_montage.py @@ -53,6 +53,28 @@ def test_shape(): assert_equal(arr_out.shape, (alpha * height, alpha * width)) +def test_rescale_intensity(): + n_images = 4 + height, width = 3, 3 + arr_in = np.arange(n_images * height * width, dtype=np.float32) + arr_in = arr_in.reshape(n_images, height, width) + + arr_out = montage2d(arr_in, rescale_intensity=True) + + gt = np.array( + [[ 0. , 0.125, 0.25 , 0. , 0.125, 0.25 ], + [ 0.375, 0.5 , 0.625, 0.375, 0.5 , 0.625], + [ 0.75 , 0.875, 1. , 0.75 , 0.875, 1. ], + [ 0. , 0.125, 0.25 , 0. , 0.125, 0.25 ], + [ 0.375, 0.5 , 0.625, 0.375, 0.5 , 0.625], + [ 0.75 , 0.875, 1. , 0.75 , 0.875, 1. ]] + ) + + assert_equal(arr_out.min(), 0.0) + assert_equal(arr_out.max(), 1.0) + assert_array_equal(arr_out, gt) + + @raises(AssertionError) def test_error_ndim(): arr_error = np.random.randn(1, 2, 3, 4) From 5e9ad08b8031752473d182a7b8f9858ee92f2e52 Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Tue, 14 Feb 2012 19:45:38 -0500 Subject: [PATCH 16/20] DOC: update view_as_blocks and gabors_from_lena example using comments from @stevanv (PR #138) --- doc/examples/plot_gabors_from_lena.py | 39 +++++++++++++-------------- doc/examples/plot_view_as_blocks.py | 27 ++++++++++--------- 2 files changed, 33 insertions(+), 33 deletions(-) diff --git a/doc/examples/plot_gabors_from_lena.py b/doc/examples/plot_gabors_from_lena.py index 0f483207..187289b7 100644 --- a/doc/examples/plot_gabors_from_lena.py +++ b/doc/examples/plot_gabors_from_lena.py @@ -7,20 +7,20 @@ Gabors / Primary Visual Cortex "Simple Cells" from Lena How to build a (bio-plausible) "sparse" dictionary (or 'codebook', or 'filterbank') for e.g. image classification without any fancy math and -with just standard python scientific librairies? +with just standard python scientific libraries? Please find below a short answer ;-) This simple example shows how to get Gabor-like filters [1]_ using just the famous Lena image. Gabor filters are good approximations of the "Simple Cells" [2]_ receptive fields [3]_ found in the mammalian primary -visual cortex (V1) (for details, see e.g. the Nobel-prize winning work of Hubel -& Wiesel done in the 60s). +visual cortex (V1) (for details, see e.g. the Nobel-prize winning work +of Hubel & Wiesel done in the 60s [4]_ [5]_). -Here we use McQueen's 'kmeans' algorithm [4]_, as a simple bio-plausible -hebbian-like learning rule and we apply it (a) to patches of the -original Lena image (retinal projection), and (b) to patches of an -LGN-like [5]_ Lena image using a simple difference of gaussians (DoG) +Here we use McQueen's 'kmeans' algorithm [6]_, as a simple biologically +plausible hebbian-like learning rule and we apply it (a) to patches of +the original Lena image (retinal projection), and (b) to patches of an +LGN-like [7]_ Lena image using a simple difference of gaussians (DoG) approximation. Enjoy ;-) And keep in mind that getting Gabors on natural image patches @@ -31,39 +31,36 @@ is not rocket science. .. [3] http://en.wikipedia.org/wiki/Receptive_field .. [4] http://en.wikipedia.org/wiki/K-means_clustering .. [5] http://en.wikipedia.org/wiki/Lateral_geniculate_nucleus - -References ----------- -D. H. Hubel and T. N. Wiesel Receptive Fields of Single Neurones in the -Cat's Striate Cortex J. Physiol. pp. 574-591 (148) 1959 - -D. H. Hubel and T. N. Wiesel Receptive Fields, Binocular Interaction and -Functional Architecture in the Cat's Visual Cortex J. Physiol. 160 pp. -106-154 1962 +.. [6] D. H. Hubel and T. N. Wiesel Receptive Fields of Single Neurones +in the Cat's Striate Cortex J. Physiol. pp. 574-591 (148) 1959 +.. [7] D. H. Hubel and T. N. Wiesel Receptive Fields, Binocular +Interaction and Functional Architecture in the Cat's Visual Cortex J. +Physiol. 160 pp. 106-154 1962 """ import numpy as np -from scipy import misc from scipy.cluster.vq import kmeans2 +from scipy import ndimage as ndi import matplotlib.pyplot as plt +from skimage import data +from skimage import color from skimage.util.shape import view_as_windows from skimage.util.montage import montage2d -from scipy import ndimage as ndi np.random.seed(42) patch_shape = 8, 8 n_filters = 49 -lena = misc.lena() / 255. +lena = color.rgb2gray(data.lena()) / 255. # -- filterbank1 on original Lena patches1 = view_as_windows(lena, patch_shape) patches1 = patches1.reshape(-1, patch_shape[0] * patch_shape[1])[::8] fb1, _ = kmeans2(patches1, n_filters, minit='points') fb1 = fb1.reshape((-1,) + patch_shape) -fb1_montage = montage2d(fb1) +fb1_montage = montage2d(fb1, rescale_intensity=True) # -- filterbank2 LGN-like Lena lena_dog = ndi.gaussian_filter(lena, .5) - ndi.gaussian_filter(lena, 1) @@ -71,7 +68,7 @@ patches2 = view_as_windows(lena_dog, patch_shape) patches2 = patches2.reshape(-1, patch_shape[0] * patch_shape[1])[::8] fb2, _ = kmeans2(patches2, n_filters, minit='points') fb2 = fb2.reshape((-1,) + patch_shape) -fb2_montage = montage2d(fb2) +fb2_montage = montage2d(fb2, rescale_intensity=True) # -- plt.figure(figsize=(9, 3)) diff --git a/doc/examples/plot_view_as_blocks.py b/doc/examples/plot_view_as_blocks.py index 305192dc..3815d3ec 100644 --- a/doc/examples/plot_view_as_blocks.py +++ b/doc/examples/plot_view_as_blocks.py @@ -3,25 +3,28 @@ Block views on images/arrays ============================ -This example illustrates the use of `view_as_blocks` from `skimage.util.shape`. -Block views can be incredibly useful when one wants to perform local operations -on non-overlapping image patches. +This example illustrates the use of `view_as_blocks` from +`skimage.util.shape`. Block views can be incredibly useful when one +wants to perform local operations on non-overlapping image patches. -We use `lena` from `scipy.misc` and virtually 'slice' it into square blocks. -Then, on each block, we either pool the mean, the max or the median value of -that block. The results are displayed altogether, along with a 'classic' -`bicubic` rescaling of the original `lena` image. +We use `lena` from `skimage.data` and virtually 'slice' it into square +blocks. Then, on each block, we either pool the mean, the max or the +median value of that block. The results are displayed altogether, along +with a 'classic' `bicubic` rescaling of the original `lena` image. """ import numpy as np -from scipy.misc import lena, imresize +from scipy.misc import imresize from matplotlib import pyplot as plt import matplotlib.cm as cm + +from skimage import data +from skimage import color from skimage.util.shape import view_as_blocks -# -- get `lena` from scipy in grayscale -l = lena() +# -- get `lena` from skimage.data in grayscale +l = color.rgb2gray(data.lena()) / 255. # -- size of blocks block_shape = (4, 4) @@ -43,7 +46,7 @@ median_view = np.median(flatten_view, axis=2) plt.figure(figsize=(10, 10)) plt.subplot(221) -plt.title("Original rescaled\n in bicubic mode"); +plt.title("Original rescaled\n in bicubic mode") l_resized = imresize(l, view.shape[:2], interp='bicubic') plt.imshow(l_resized, cmap=cm.Greys_r) @@ -52,7 +55,7 @@ plt.title("Block view with\n local mean pooling") plt.imshow(mean_view, cmap=cm.Greys_r) plt.subplot(223) -plt.title("Block view with\n local max pooling"); +plt.title("Block view with\n local max pooling") plt.imshow(max_view, cmap=cm.Greys_r) plt.subplot(224) From e77d7a7e443b6b268e690cbb64051148a4b3ca0f Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Tue, 14 Feb 2012 19:56:42 -0500 Subject: [PATCH 17/20] DOC: don't divide by 255 when it's not required ;-) --- doc/examples/plot_gabors_from_lena.py | 2 +- doc/examples/plot_view_as_blocks.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/examples/plot_gabors_from_lena.py b/doc/examples/plot_gabors_from_lena.py index 187289b7..36872259 100644 --- a/doc/examples/plot_gabors_from_lena.py +++ b/doc/examples/plot_gabors_from_lena.py @@ -53,7 +53,7 @@ np.random.seed(42) patch_shape = 8, 8 n_filters = 49 -lena = color.rgb2gray(data.lena()) / 255. +lena = color.rgb2gray(data.lena()) # -- filterbank1 on original Lena patches1 = view_as_windows(lena, patch_shape) diff --git a/doc/examples/plot_view_as_blocks.py b/doc/examples/plot_view_as_blocks.py index 3815d3ec..41388a89 100644 --- a/doc/examples/plot_view_as_blocks.py +++ b/doc/examples/plot_view_as_blocks.py @@ -24,7 +24,7 @@ from skimage.util.shape import view_as_blocks # -- get `lena` from skimage.data in grayscale -l = color.rgb2gray(data.lena()) / 255. +l = color.rgb2gray(data.lena()) # -- size of blocks block_shape = (4, 4) From ac54f6969e13e4492db2471cb9465a5c6d82e92b Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Tue, 14 Feb 2012 20:01:10 -0500 Subject: [PATCH 18/20] DOC: remove dependency on PIL, see PR #138 --- doc/examples/plot_view_as_blocks.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/doc/examples/plot_view_as_blocks.py b/doc/examples/plot_view_as_blocks.py index 41388a89..876e5ab9 100644 --- a/doc/examples/plot_view_as_blocks.py +++ b/doc/examples/plot_view_as_blocks.py @@ -10,11 +10,12 @@ wants to perform local operations on non-overlapping image patches. We use `lena` from `skimage.data` and virtually 'slice' it into square blocks. Then, on each block, we either pool the mean, the max or the median value of that block. The results are displayed altogether, along -with a 'classic' `bicubic` rescaling of the original `lena` image. +with a spline interpolation of order 3 rescaling of the original `lena` +image. """ import numpy as np -from scipy.misc import imresize +from scipy import ndimage as ndi from matplotlib import pyplot as plt import matplotlib.cm as cm @@ -46,8 +47,8 @@ median_view = np.median(flatten_view, axis=2) plt.figure(figsize=(10, 10)) plt.subplot(221) -plt.title("Original rescaled\n in bicubic mode") -l_resized = imresize(l, view.shape[:2], interp='bicubic') +plt.title("Original rescaled with\n spline interpolation (order=3)") +l_resized = ndi.zoom(l, 2, order=3) plt.imshow(l_resized, cmap=cm.Greys_r) plt.subplot(222) From c61fe51d2f5edd44f4a7be3b81d18840dfca844d Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Tue, 14 Feb 2012 20:08:25 -0500 Subject: [PATCH 19/20] DOC: remove nastiness, replace with appropriate doctest directive --- skimage/util/montage.py | 4 +--- skimage/util/shape.py | 15 +++------------ 2 files changed, 4 insertions(+), 15 deletions(-) diff --git a/skimage/util/montage.py b/skimage/util/montage.py index a879a0a2..7d9fe5ee 100644 --- a/skimage/util/montage.py +++ b/skimage/util/montage.py @@ -50,13 +50,11 @@ def montage2d(arr_in, fill='mean', rescale_intensity=False): >>> import numpy as np >>> from skimage.util.montage import montage2d >>> arr_in = np.arange(3 * 2 * 2).reshape(3, 2, 2) - >>> print arr_in + >>> print arr_in # doctest: +NORMALIZE_WHITESPACE [[[ 0 1] [ 2 3]] - [[ 4 5] [ 6 7]] - [[ 8 9] [10 11]]] >>> arr_out = montage2d(arr_in) diff --git a/skimage/util/shape.py b/skimage/util/shape.py index 7d796405..fc562ea9 100644 --- a/skimage/util/shape.py +++ b/skimage/util/shape.py @@ -37,22 +37,19 @@ def view_as_blocks(arr_in, block_shape): >>> B[1, 0, 1, 1] 13 >>> A = np.arange(4*4*6).reshape(4,4,6) - >>> A + >>> A # doctest: +NORMALIZE_WHITESPACE 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], @@ -60,11 +57,9 @@ def view_as_blocks(arr_in, block_shape): >>> B = view_as_blocks(A, block_shape=(1, 2, 2)) >>> B.shape (4, 2, 3, 1, 2, 2) - >>> B[2:, 0, 2] + >>> B[2:, 0, 2] # doctest: +NORMALIZE_WHITESPACE array([[[[52, 53], [58, 59]]], - - [[[76, 77], [82, 83]]]]) """ @@ -166,23 +161,19 @@ def view_as_windows(arr_in, window_shape): >>> B = view_as_windows(A, window_shape) >>> B.shape (2, 2, 4, 3) - >>> B + >>> B # doctest: +NORMALIZE_WHITESPACE 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], From 0a9522d4ea011698c6a94c369271ec61ddc35e24 Mon Sep 17 00:00:00 2001 From: Nicolas Pinto Date: Tue, 14 Feb 2012 21:56:48 -0500 Subject: [PATCH 20/20] FIX: montage2d should not use rescale_intensity's out_range kwarg --- skimage/util/montage.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/skimage/util/montage.py b/skimage/util/montage.py index 7d9fe5ee..fe3c68a2 100644 --- a/skimage/util/montage.py +++ b/skimage/util/montage.py @@ -75,8 +75,7 @@ def montage2d(arr_in, fill='mean', rescale_intensity=False): # -- rescale intensity if necessary if rescale_intensity: for i in xrange(n_images): - arr_in[i] = exposure.rescale_intensity( - arr_in[i], out_range=(0.0, 1.0)) + arr_in[i] = exposure.rescale_intensity(arr_in[i]) # -- determine alpha alpha = int(np.ceil(np.sqrt(n_images)))