diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 8f479945..5082c508 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -2,7 +2,9 @@ Project coordination - Nicolas Pinto - Colour spaces and filters + 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 @@ -93,3 +95,6 @@ - Gaël Varoquaux Harris corner detector + +- Nicolas Poilvert + Shape views: ``util.shape.view_as_windows`` and ``util.shape.view_as_blocks`` diff --git a/doc/examples/plot_gabors_from_lena.py b/doc/examples/plot_gabors_from_lena.py new file mode 100644 index 00000000..36872259 --- /dev/null +++ b/doc/examples/plot_gabors_from_lena.py @@ -0,0 +1,97 @@ +""" +======================================================= +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 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 [4]_ [5]_). + +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 +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 +.. [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.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 + +np.random.seed(42) + +patch_shape = 8, 8 +n_filters = 49 + +lena = color.rgb2gray(data.lena()) + +# -- 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, rescale_intensity=True) + +# -- 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, rescale_intensity=True) + +# -- +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() diff --git a/doc/examples/plot_view_as_blocks.py b/doc/examples/plot_view_as_blocks.py new file mode 100644 index 00000000..876e5ab9 --- /dev/null +++ b/doc/examples/plot_view_as_blocks.py @@ -0,0 +1,67 @@ +""" +============================ +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 `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 spline interpolation of order 3 rescaling of the original `lena` +image. +""" + +import numpy as np +from scipy import ndimage as ndi +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 skimage.data in grayscale +l = color.rgb2gray(data.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 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) +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() diff --git a/skimage/util/__init__.py b/skimage/util/__init__.py index 6ef65e20..980e3880 100644 --- a/skimage/util/__init__.py +++ b/skimage/util/__init__.py @@ -1,2 +1 @@ from .dtype import * - diff --git a/skimage/util/montage.py b/skimage/util/montage.py new file mode 100644 index 00000000..fe3c68a2 --- /dev/null +++ b/skimage/util/montage.py @@ -0,0 +1,96 @@ +__all__ = ['montage2d'] + +import numpy as np +from .. import exposure + +EPSILON = 1e-6 + + +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. + + 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', 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(). + + rescale_intensity: bool, optional + Whether to rescale the intensity of each image to [0, 1]. + + 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 # doctest: +NORMALIZE_WHITESPACE + [[[ 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 + + # -- rescale intensity if necessary + if rescale_intensity: + for i in xrange(n_images): + arr_in[i] = exposure.rescale_intensity(arr_in[i]) + + # -- 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/shape.py b/skimage/util/shape.py new file mode 100644 index 00000000..fc562ea9 --- /dev/null +++ b/skimage/util/shape.py @@ -0,0 +1,208 @@ +__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-dimensional 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 # 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], + [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] # doctest: +NORMALIZE_WHITESPACE + 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 # 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], + [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_montage.py b/skimage/util/tests/test_montage.py new file mode 100644 index 00000000..47e5426d --- /dev/null +++ b/skimage/util/tests/test_montage.py @@ -0,0 +1,81 @@ +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)) + + +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) + montage2d(arr_error) diff --git a/skimage/util/tests/test_shape.py b/skimage/util/tests/test_shape.py new file mode 100644 index 00000000..3b5f8d41 --- /dev/null +++ b/skimage/util/tests/test_shape.py @@ -0,0 +1,141 @@ +import numpy as np +from nose.tools import raises +from numpy.testing import assert_equal +from skimage.util.shape import view_as_blocks, view_as_windows + + +@raises(TypeError) +def test_view_as_blocks_block_not_a_tuple(): + + A = np.arange(10) + view_as_blocks(A, [5]) + + +@raises(ValueError) +def test_view_as_blocks_negative_shape(): + + A = np.arange(10) + view_as_blocks(A, (-2,)) + + +@raises(ValueError) +def test_view_as_blocks_block_too_large(): + + A = np.arange(10) + view_as_blocks(A, (11,)) + + +@raises(ValueError) +def test_view_as_blocks_wrong_block_dimension(): + + A = np.arange(10) + view_as_blocks(A, (2, 2)) + + +@raises(ValueError) +def test_view_as_blocks_1D_array_wrong_block_shape(): + + A = np.arange(10) + view_as_blocks(A, (3,)) + + +def test_view_as_blocks_1D_array(): + + A = np.arange(10) + B = view_as_blocks(A, (5,)) + assert_equal(B, np.array([[0, 1, 2, 3, 4], + [5, 6, 7, 8, 9]])) + + +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], + [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)) + 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]]]])) + + +@raises(TypeError) +def test_view_as_windows_input_not_array(): + + A = [1, 2, 3, 4, 5] + view_as_windows(A, (2,)) + + +@raises(TypeError) +def test_view_as_windows_window_not_tuple(): + + A = np.arange(10) + view_as_windows(A, [2]) + + +@raises(ValueError) +def test_view_as_windows_wrong_window_dimension(): + + A = np.arange(10) + view_as_windows(A, (2, 2)) + + +@raises(ValueError) +def test_view_as_windows_negative_window_length(): + + A = np.arange(10) + view_as_windows(A, (-1,)) + + +@raises(ValueError) +def test_view_as_windows_window_too_large(): + + A = np.arange(10) + view_as_windows(A, (11,)) + + +def test_view_as_windows_1D(): + + A = np.arange(10) + window_shape = (3,) + B = view_as_windows(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_view_as_windows_2D(): + + 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], + [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]]]]))