wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-20 12:40:31 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #138 from poilvert/array_views

Browse Source 
ENH: Add window views and montage.
This commit is contained in:
Stefan van der Walt
2012-02-19 00:27:09 -08:00
parent e1e3d6ce2b 0a9522d4ea
commit 0c304d5acd
8 changed files with 696 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CONTRIBUTORS.txt
+6 -1
View File
@@ -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``
doc/examples/plot_gabors_from_lena.py
+97
View File
@@ -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()
doc/examples/plot_view_as_blocks.py
+67
View File
@@ -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()
skimage/util/__init__.py
-1
View File
@@ -1,2 +1 @@
from .dtype import *
skimage/util/montage.py
+96
View File
@@ -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
skimage/util/shape.py
+208
View File
@@ -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
skimage/util/tests/test_montage.py
+81
View File
@@ -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)
skimage/util/tests/test_shape.py
+141
View File
@@ -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]]]]))