From 56764f82095afb0a32408b224ef23c660d4d87b0 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 18 Dec 2014 20:02:08 +1100 Subject: [PATCH 01/36] Wrap scipy.ndimage grey morphology functions --- skimage/morphology/grey.py | 141 ++++++++++++++++++++++--------------- skimage/morphology/misc.py | 26 ++----- 2 files changed, 90 insertions(+), 77 deletions(-) diff --git a/skimage/morphology/grey.py b/skimage/morphology/grey.py index e65e756b..0dbe183d 100644 --- a/skimage/morphology/grey.py +++ b/skimage/morphology/grey.py @@ -1,17 +1,53 @@ """ Grayscale morphological operations """ -from skimage import img_as_ubyte -from .misc import default_fallback - -from . import cmorph - +import numpy as np +from scipy import ndimage as nd +from .misc import default_selem __all__ = ['erosion', 'dilation', 'opening', 'closing', 'white_tophat', 'black_tophat'] -@default_fallback +def _shift_selem(selem, shift_x, shift_y): + """Shift the binary image `selem` in the left and/or up. + + This only affects structuring elements with even number of rows or + columns. + + Parameters + ---------- + selem : 2D array, shape (M, N) + The input structuring element. + shift_x, shift_y : bool + Whether to move `selem` along each axis. + + Returns + ------- + out : 2D array, shape (M + int(shift_x), N + int(shift_y)) + The shifted structuring element. + """ + if selem.ndim > 2: + # do nothing for 3D or higher structuring elements + return selem + m, n = selem.shape + if m % 2 == 0: + extra_row = np.zeros((1, n), selem.dtype) + if shift_x: + selem = np.vstack((selem, extra_row)) + else: + selem = np.vstack((extra_row, selem)) + m += 1 + if n % 2 == 0: + extra_col = np.zeros((m, 1), selem.dtype) + if shift_y: + selem = np.hstack((selem, extra_col)) + else: + selem = np.hstack((extra_col, selem)) + return selem + + +@default_selem def erosion(image, selem=None, out=None, shift_x=False, shift_y=False): """Return greyscale morphological erosion of an image. @@ -24,7 +60,7 @@ def erosion(image, selem=None, out=None, shift_x=False, shift_y=False): image : ndarray Image array. selem : ndarray, optional - The neighborhood expressed as a 2-D array of 1's and 0's. + The neighborhood expressed as an array of 1's and 0's. If None, use cross-shaped structuring element (connectivity=1). out : ndarrays, optional The array to store the result of the morphology. If None is @@ -35,7 +71,7 @@ def erosion(image, selem=None, out=None, shift_x=False, shift_y=False): Returns ------- - eroded : uint8 array + eroded : array, same shape as `image` The result of the morphological erosion. Notes @@ -62,16 +98,15 @@ def erosion(image, selem=None, out=None, shift_x=False, shift_y=False): [0, 0, 0, 0, 0]], dtype=uint8) """ - - if image is out: - raise NotImplementedError("In-place erosion not supported!") - image = img_as_ubyte(image) - selem = img_as_ubyte(selem) - return cmorph._erode(image, selem, out=out, - shift_x=shift_x, shift_y=shift_y) + selem = np.array(selem) + selem = _shift_selem(selem, shift_x, shift_y) + if out is None: + out = np.empty_like(image) + nd.grey_erosion(image, footprint=selem, output=out) + return out -@default_fallback +@default_selem def dilation(image, selem=None, out=None, shift_x=False, shift_y=False): """Return greyscale morphological dilation of an image. @@ -123,17 +158,15 @@ def dilation(image, selem=None, out=None, shift_x=False, shift_y=False): [0, 0, 0, 0, 0]], dtype=uint8) """ - - if image is out: - raise NotImplementedError("In-place dilation not supported!") - - image = img_as_ubyte(image) - selem = img_as_ubyte(selem) - return cmorph._dilate(image, selem, out=out, - shift_x=shift_x, shift_y=shift_y) + selem = np.array(selem) + selem = _shift_selem(selem, shift_x, shift_y) + if out is None: + out = np.empty_like(image) + nd.grey_dilation(image, footprint=selem, output=out) + return out -@default_fallback +@default_selem def opening(image, selem=None, out=None): """Return greyscale morphological opening of an image. @@ -147,7 +180,7 @@ def opening(image, selem=None, out=None): image : ndarray Image array. selem : ndarray, optional - The neighborhood expressed as a 2-D array of 1's and 0's. + The neighborhood expressed as an array of 1's and 0's. If None, use cross-shaped structuring element (connectivity=1). out : ndarray, optional The array to store the result of the morphology. If None @@ -155,7 +188,7 @@ def opening(image, selem=None, out=None): Returns ------- - opening : uint8 array + opening : array The result of the morphological opening. Examples @@ -176,17 +209,13 @@ def opening(image, selem=None, out=None): [0, 0, 0, 0, 0]], dtype=uint8) """ - - h, w = selem.shape - shift_x = True if (w % 2) == 0 else False - shift_y = True if (h % 2) == 0 else False - eroded = erosion(image, selem) - out = dilation(eroded, selem, out=out, shift_x=shift_x, shift_y=shift_y) + # note: shift_x, shift_y do nothing if selem side length is odd + out = dilation(eroded, selem, out=out, shift_x=True, shift_y=True) return out -@default_fallback +@default_selem def closing(image, selem=None, out=None): """Return greyscale morphological closing of an image. @@ -200,7 +229,7 @@ def closing(image, selem=None, out=None): image : ndarray Image array. selem : ndarray, optional - The neighborhood expressed as a 2-D array of 1's and 0's. + The neighborhood expressed as an array of 1's and 0's. If None, use cross-shaped structuring element (connectivity=1). out : ndarray, optional The array to store the result of the morphology. If None, @@ -208,7 +237,7 @@ def closing(image, selem=None, out=None): Returns ------- - closing : uint8 array + closing : array The result of the morphological closing. Examples @@ -229,17 +258,13 @@ def closing(image, selem=None, out=None): [0, 0, 0, 0, 0]], dtype=uint8) """ - - h, w = selem.shape - shift_x = True if (w % 2) == 0 else False - shift_y = True if (h % 2) == 0 else False - dilated = dilation(image, selem) - out = erosion(dilated, selem, out=out, shift_x=shift_x, shift_y=shift_y) + # note: shift_x, shift_y do nothing if selem side length is odd + out = erosion(dilated, selem, out=out, shift_x=True, shift_y=True) return out -@default_fallback +@default_selem def white_tophat(image, selem=None, out=None): """Return white top hat of an image. @@ -252,7 +277,7 @@ def white_tophat(image, selem=None, out=None): image : ndarray Image array. selem : ndarray, optional - The neighborhood expressed as a 2-D array of 1's and 0's. + The neighborhood expressed as an array of 1's and 0's. If None, use cross-shaped structuring element (connectivity=1). out : ndarray, optional The array to store the result of the morphology. If None @@ -260,7 +285,7 @@ def white_tophat(image, selem=None, out=None): Returns ------- - opening : uint8 array + out : array The result of the morphological white top hat. Examples @@ -281,16 +306,18 @@ def white_tophat(image, selem=None, out=None): [0, 0, 0, 0, 0]], dtype=uint8) """ - - if image is out: - raise NotImplementedError("Cannot perform white top hat in place.") - - out = opening(image, selem, out=out) - out = image - out + selem = np.array(selem) + if out is image: + opened = opening(image, selem) + out -= opened + return out + elif out is None: + out = np.empty_like(image) + out = nd.white_tophat(image, footprint=selem, output=out) return out -@default_fallback +@default_selem def black_tophat(image, selem=None, out=None): """Return black top hat of an image. @@ -333,10 +360,10 @@ def black_tophat(image, selem=None, out=None): [0, 0, 0, 0, 0]], dtype=uint8) """ - - if image is out: - raise NotImplementedError("Cannot perform white top hat in place.") - + if out is image: + original = image.copy() + else: + original = image out = closing(image, selem, out=out) - out = out - image + out -= original return out diff --git a/skimage/morphology/misc.py b/skimage/morphology/misc.py index 55b279a2..663e1e02 100644 --- a/skimage/morphology/misc.py +++ b/skimage/morphology/misc.py @@ -15,11 +15,8 @@ funcs = ('binary_erosion', 'binary_dilation', 'binary_opening', skimage2ndimage.update(dict((x, x) for x in funcs)) -def default_fallback(func): - """Decorator to fall back on ndimage for images with more than 2 dimensions - - Decorator also provides a default structuring element, `selem`, with the - appropriate dimensionality if none is specified. +def default_selem(func): + """Decorator to add a default structuring element to morphology functions. Parameters ---------- @@ -30,25 +27,14 @@ def default_fallback(func): Returns ------- func_out : function - If the image dimensionality is greater than 2D, the ndimage - function is returned, otherwise skimage function is used. + The function, using a default structuring element of same dimension + as the input image with connectivity 1. """ @functools.wraps(func) - def func_out(image, selem=None, out=None, **kwargs): - # Default structure element + def func_out(image, selem=None, *args, **kwargs): if selem is None: selem = _default_selem(image.ndim) - - # If image has more than 2 dimensions, use scipy.ndimage - if image.ndim > 2: - function = getattr(nd, skimage2ndimage[func.__name__]) - try: - return function(image, footprint=selem, output=out, **kwargs) - except TypeError: - # nd.binary_* take structure instead of footprint - return function(image, structure=selem, output=out, **kwargs) - else: - return func(image, selem=selem, out=out, **kwargs) + return func(image, selem=selem, *args, **kwargs) return func_out From 8ffed5e424f48adea60b8bd1bb16c639c522e9a8 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 18 Dec 2014 21:09:42 +1100 Subject: [PATCH 02/36] Wrap scipy.ndimage binary morphology functions --- skimage/morphology/binary.py | 52 ++++++++++-------------------------- 1 file changed, 14 insertions(+), 38 deletions(-) diff --git a/skimage/morphology/binary.py b/skimage/morphology/binary.py index 4244738f..cc60d047 100644 --- a/skimage/morphology/binary.py +++ b/skimage/morphology/binary.py @@ -2,15 +2,13 @@ Binary morphological operations """ import numpy as np -from scipy import ndimage -from .misc import default_fallback +from scipy import ndimage as nd +from .misc import default_selem -# Our functions only work in 2D, so for 3D or higher input we should fall back -# on `scipy.ndimage`. Additionally, we want to use a cross-shaped structuring -# element of the appropriate dimension for each of these functions. -# The `default_callback` provides all these. -@default_fallback +# The default_selem decorator provides a diamond structuring element as default +# with the appropriate dimension for the input `image`. +@default_selem def binary_erosion(image, selem=None, out=None): """Return fast binary morphological erosion of an image. @@ -38,24 +36,13 @@ def binary_erosion(image, selem=None, out=None): The result of the morphological erosion with values in ``[0, 1]``. """ - - selem = (selem != 0) - selem_sum = np.sum(selem) - - if selem_sum <= 255: - conv = np.empty_like(image, dtype=np.uint8) - else: - conv = np.empty_like(image, dtype=np.uint) - - binary = (image > 0).view(np.uint8) - ndimage.convolve(binary, selem, mode='constant', cval=1, output=conv) - if out is None: - out = np.empty_like(conv, dtype=np.bool) - return np.equal(conv, selem_sum, out=out) + out = np.empty(image.shape, dtype=np.bool) + nd.binary_erosion(image, structure=selem, output=out) + return out -@default_fallback +@default_selem def binary_dilation(image, selem=None, out=None): """Return fast binary morphological dilation of an image. @@ -84,23 +71,13 @@ def binary_dilation(image, selem=None, out=None): The result of the morphological dilation with values in ``[0, 1]``. """ - - selem = (selem != 0) - - if np.sum(selem) <= 255: - conv = np.empty_like(image, dtype=np.uint8) - else: - conv = np.empty_like(image, dtype=np.uint) - - binary = (image > 0).view(np.uint8) - ndimage.convolve(binary, selem, mode='constant', cval=0, output=conv) - if out is None: - out = np.empty_like(conv, dtype=np.bool) - return np.not_equal(conv, 0, out=out) + out = np.empty(image.shape, dtype=np.bool) + nd.binary_dilation(image, structure=selem, output=out) + return out -@default_fallback +@default_selem def binary_opening(image, selem=None, out=None): """Return fast binary morphological opening of an image. @@ -134,7 +111,7 @@ def binary_opening(image, selem=None, out=None): return out -@default_fallback +@default_selem def binary_closing(image, selem=None, out=None): """Return fast binary morphological closing of an image. @@ -163,7 +140,6 @@ def binary_closing(image, selem=None, out=None): The result of the morphological closing. """ - dilated = binary_dilation(image, selem) out = binary_erosion(dilated, selem, out=out) return out From 1fa6904c04715a4be0c37a30d5552170400b293c Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 18 Dec 2014 22:33:17 +1100 Subject: [PATCH 03/36] Workaround for strange ndimage selem inversion --- skimage/morphology/grey.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/skimage/morphology/grey.py b/skimage/morphology/grey.py index 0dbe183d..c4ff5ffd 100644 --- a/skimage/morphology/grey.py +++ b/skimage/morphology/grey.py @@ -47,6 +47,38 @@ def _shift_selem(selem, shift_x, shift_y): return selem +def _invert_selem(selem): + """Change the order of the values in `selem`. + + This is a patch for the *weird* footprint inversion in + `nd.grey_morphology` [1]. + + Parameters + ---------- + selem : array + The input structuring element. + + Returns + ------- + inverted : array + The structuring element, in opposite order. + + Examples + -------- + >>> selem = np.array([[0, 0, 0], [0, 1, 1], [0, 1, 1]], np.uint8) + >>> _invert_selem(selem) + array([[1, 1, 0], + [1, 1, 0], + [0, 0, 0]], dtype=uint8) + + References + ---------- + [1] https://github.com/scipy/scipy/blob/ec20ababa400e39ac3ffc9148c01ef86d5349332/scipy/ndimage/morphology.py#L1285 + """ + inverted = selem[(slice(None, None, -1),) * selem.ndim] + return inverted + + @default_selem def erosion(image, selem=None, out=None, shift_x=False, shift_y=False): """Return greyscale morphological erosion of an image. @@ -160,6 +192,10 @@ def dilation(image, selem=None, out=None, shift_x=False, shift_y=False): """ selem = np.array(selem) selem = _shift_selem(selem, shift_x, shift_y) + # invert the structuring element to patch the same thing happening inside + # ndimage.grey_dilation + # https://github.com/scipy/scipy/blob/ec20ababa400e39ac3ffc9148c01ef86d5349332/scipy/ndimage/morphology.py#L1285 + selem = _invert_selem(selem) if out is None: out = np.empty_like(image) nd.grey_dilation(image, footprint=selem, output=out) From 04536cc7dfe392a458436efe80cfd0010faad537 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 5 Jan 2015 18:34:46 +1100 Subject: [PATCH 04/36] Update morphology dtypes tests Previously, we were testing that any dtype would get converted to uint8 and then correctly processed. Now, since we are using ndimage, we are directly processing all dtypes. I've updated the tests accordingly. --- skimage/morphology/tests/test_grey.py | 78 +++++++++++++-------------- 1 file changed, 39 insertions(+), 39 deletions(-) diff --git a/skimage/morphology/tests/test_grey.py b/skimage/morphology/tests/test_grey.py index 911e0c71..576404f2 100644 --- a/skimage/morphology/tests/test_grey.py +++ b/skimage/morphology/tests/test_grey.py @@ -5,7 +5,7 @@ from numpy import testing from scipy import ndimage import skimage -from skimage import data_dir +from skimage import data_dir, img_as_uint from skimage.morphology import grey, selem from skimage._shared._warnings import expected_warnings @@ -208,51 +208,51 @@ def test_2d_ndimage_equivalence(): testing.assert_array_equal(opened, ndimage_opened) testing.assert_array_equal(closed, ndimage_closed) -class TestDTypes(): +# float test images +im = np.array([[ 0.55, 0.72, 0.6 , 0.54, 0.42], + [ 0.65, 0.44, 0.89, 0.96, 0.38], + [ 0.79, 0.53, 0.57, 0.93, 0.07], + [ 0.09, 0.02, 0.83, 0.78, 0.87], + [ 0.98, 0.8 , 0.46, 0.78, 0.12]]) - def setUp(self): - k = 5 - arrname = '%03i' % k +eroded = np.array([[ 0.55, 0.44, 0.54, 0.42, 0.38], + [ 0.44, 0.44, 0.44, 0.38, 0.07], + [ 0.09, 0.02, 0.53, 0.07, 0.07], + [ 0.02, 0.02, 0.02, 0.78, 0.07], + [ 0.09, 0.02, 0.46, 0.12, 0.12]]) - self.disk = selem.disk(k) +dilated = np.array([[ 0.72, 0.72, 0.89, 0.96, 0.54], + [ 0.79, 0.89, 0.96, 0.96, 0.96], + [ 0.79, 0.79, 0.93, 0.96, 0.93], + [ 0.98, 0.83, 0.83, 0.93, 0.87], + [ 0.98, 0.98, 0.83, 0.78, 0.87]]) - fname_opening = os.path.join(data_dir, "disk-open-matlab-output.npz") - self.expected_opening = np.load(fname_opening)[arrname] +opened = np.array([[ 0.55, 0.55, 0.54, 0.54, 0.42], + [ 0.55, 0.44, 0.54, 0.44, 0.38], + [ 0.44, 0.53, 0.53, 0.78, 0.07], + [ 0.09, 0.02, 0.78, 0.78, 0.78], + [ 0.09, 0.46, 0.46, 0.78, 0.12]]) - fname_closing = os.path.join(data_dir, "disk-close-matlab-output.npz") - self.expected_closing = np.load(fname_closing)[arrname] +closed = np.array([[ 0.72, 0.72, 0.72, 0.54, 0.54], + [ 0.72, 0.72, 0.89, 0.96, 0.54], + [ 0.79, 0.79, 0.79, 0.93, 0.87], + [ 0.79, 0.79, 0.83, 0.78, 0.87], + [ 0.98, 0.83, 0.78, 0.78, 0.78]]) - def _test_image(self, image): - with expected_warnings(['precision loss']): - result_opening = grey.opening(image, self.disk) - testing.assert_equal(result_opening, self.expected_opening) - - with expected_warnings(['precision loss']): - result_closing = grey.closing(image, self.disk) - testing.assert_equal(result_closing, self.expected_closing) - - def test_float(self): - image = skimage.img_as_float(lena) - self._test_image(image) - - @testing.decorators.skipif(True) - def test_int(self): - image = skimage.img_as_int(lena) - self._test_image(image) - - def test_uint(self): - image = skimage.img_as_uint(lena) - self._test_image(image) +def test_float(): + np.testing.assert_allclose(grey.erosion(im), eroded) + np.testing.assert_allclose(grey.dilation(im), dilated) + np.testing.assert_allclose(grey.opening(im), opened) + np.testing.assert_allclose(grey.closing(im), closed) -def test_inplace(): - selem = np.ones((3, 3)) - image = np.zeros((5, 5)) - out = image - - for f in (grey.erosion, grey.dilation, - grey.white_tophat, grey.black_tophat): - testing.assert_raises(NotImplementedError, f, image, selem, out=out) +def test_uint16(): + im16, eroded16, dilated16, opened16, closed16 = ( + map(img_as_uint, [im, eroded, dilated, opened, closed])) + np.testing.assert_allclose(grey.erosion(im16), eroded16) + np.testing.assert_allclose(grey.dilation(im16), dilated16) + np.testing.assert_allclose(grey.opening(im16), opened16) + np.testing.assert_allclose(grey.closing(im16), closed16) def test_discontiguous_out_array(): From 97770e9fa5bd282550b4ef4d0d3bd6a009cfb5f9 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 5 Jan 2015 18:48:47 +1100 Subject: [PATCH 05/36] Minor update to binary morphology docstrings --- skimage/morphology/binary.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/skimage/morphology/binary.py b/skimage/morphology/binary.py index cc60d047..b12b2948 100644 --- a/skimage/morphology/binary.py +++ b/skimage/morphology/binary.py @@ -33,7 +33,8 @@ def binary_erosion(image, selem=None, out=None): Returns ------- eroded : ndarray of bool or uint - The result of the morphological erosion with values in ``[0, 1]``. + The result of the morphological erosion taking values in + ``[False, True]``. """ if out is None: @@ -68,8 +69,8 @@ def binary_dilation(image, selem=None, out=None): Returns ------- dilated : ndarray of bool or uint - The result of the morphological dilation with values in ``[0, 1]``. - + The result of the morphological dilation with values in + ``[False, True]``. """ if out is None: out = np.empty(image.shape, dtype=np.bool) From fb7ed1d2532b55961b24224260c50988e61b4da6 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Tue, 6 Jan 2015 11:27:27 +1100 Subject: [PATCH 06/36] Add crop function, inverse to pad --- skimage/util/arraypad.py | 42 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 41 insertions(+), 1 deletion(-) diff --git a/skimage/util/arraypad.py b/skimage/util/arraypad.py index d616102d..ee0a3c5d 100644 --- a/skimage/util/arraypad.py +++ b/skimage/util/arraypad.py @@ -8,7 +8,7 @@ from __future__ import division, absolute_import, print_function import numpy as np -__all__ = ['pad'] +__all__ = ['pad', 'crop'] ############################################################################### @@ -1496,3 +1496,43 @@ def pad(array, pad_width, mode=None, **kwargs): newmat = _pad_wrap(newmat, (pad_before, pad_after), axis) return newmat + + +def crop(ar, crop_width, copy=False, order='K'): + """Crop array `ar` by `crop_width` along each dimension. + + Parameters + ---------- + ar : array-like of rank N + Input array. + crop_width : {sequence, int} + Number of values to remove from the edges of each axis. + ((before_1, after_1), ... (before_N, after_N)) specifies unique + crop widths at the start and end of each axis. + ((before, after),) specifies the same start and end crops for + every axis. + (int,) or int is a shortcut for before = after = int for all + axes. + copy : bool, optional + Ensure that the returned array is contiguous. Normally, a crop + operation will return a discontiguous view of the underlying + input array. Passing `copy=True` will result in a contiguous + copy. + order : {'C', 'F', 'A', 'K'}, optional + If `copy==True`, control the memory layout of the copy. See + `np.copy`. + + Returns + ------- + cropped : array + The cropped array. If `copy=False` (default), this is a sliced + view of the input array. + """ + ar = np.array(ar) + crops = _validate_lengths(ar, crop_width) + slices = tuple([slice(a, -b) for a, b in crops]) + if copy: + cropped = np.copy(ar[slices], order=order) + else: + cropped = ar[slices] + return cropped From eaa2aec3cc3423c60cdad604b054af1a58b7b810 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Tue, 6 Jan 2015 17:58:18 +1100 Subject: [PATCH 07/36] Add test functions for crop --- skimage/util/tests/test_arraypad.py | 42 +++++++++++++++++++++++++++-- 1 file changed, 40 insertions(+), 2 deletions(-) diff --git a/skimage/util/tests/test_arraypad.py b/skimage/util/tests/test_arraypad.py index eb8554b1..df4f0b5e 100644 --- a/skimage/util/tests/test_arraypad.py +++ b/skimage/util/tests/test_arraypad.py @@ -5,8 +5,8 @@ from __future__ import division, absolute_import, print_function import numpy as np from numpy.testing import (assert_array_equal, assert_raises, assert_allclose, - TestCase) -from skimage.util import pad + assert_equal, TestCase) +from skimage.util import pad, crop class TestConditionalShortcuts(TestCase): @@ -1043,5 +1043,43 @@ class TypeError1(TestCase): **kwargs) +def test_multi_crop(): + arr = np.arange(45).reshape(9, 5) + out = crop(arr, (1, 2), (2, 1)) + assert_array_equal(out[0], [7, 8]) + assert_array_equal(out[-1], [32, 33]) + assert_equal(out.shape, (6, 2)) + + +def test_pair_crop(): + arr = np.arange(45).reshape(9, 5) + out = crop(arr, (1, 2)) + assert_array_equal(out[0], [6, 7]) + assert_array_equal(out[-1], [31, 32]) + assert_equal(out.shape, (6, 2)) + + +def test_int_crop(): + arr = np.arange(45).reshape(9, 5) + out = crop(arr, 1) + assert_array_equal(out[0], [6, 7, 8]) + assert_array_equal(out[-1], [36, 37, 38]) + assert_equal(out.shape, (7, 3)) + + +def test_copy_crop(): + arr = np.arange(45).reshape(9, 5) + out0 = crop(arr, 1, copy=True) + assert out0.c_contiguous + out0[0, 0] = 100 + assert not np.any(arr == 100) + assert not np.may_share_memory(arr, out0) + + out1 = crop(arr, 1) + out1[0, 0] = 100 + assert arr[1, 1] == 100 + assert np.may_share_memory(arr, out1) + + if __name__ == "__main__": np.testing.run_module_suite() From d7138f1a515f7b721c1472d979a49315c0a15006 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 7 Jan 2015 11:24:45 +1100 Subject: [PATCH 08/36] Fix incorrect crop tests --- skimage/util/tests/test_arraypad.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/util/tests/test_arraypad.py b/skimage/util/tests/test_arraypad.py index df4f0b5e..4d990d0a 100644 --- a/skimage/util/tests/test_arraypad.py +++ b/skimage/util/tests/test_arraypad.py @@ -1045,7 +1045,7 @@ class TypeError1(TestCase): def test_multi_crop(): arr = np.arange(45).reshape(9, 5) - out = crop(arr, (1, 2), (2, 1)) + out = crop(arr, ((1, 2), (2, 1))) assert_array_equal(out[0], [7, 8]) assert_array_equal(out[-1], [32, 33]) assert_equal(out.shape, (6, 2)) @@ -1070,7 +1070,7 @@ def test_int_crop(): def test_copy_crop(): arr = np.arange(45).reshape(9, 5) out0 = crop(arr, 1, copy=True) - assert out0.c_contiguous + assert out0.flags.c_contiguous out0[0, 0] = 100 assert not np.any(arr == 100) assert not np.may_share_memory(arr, out0) From c96fdb9338106546326fcd71cc0c067b6506038d Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 7 Jan 2015 11:25:30 +1100 Subject: [PATCH 09/36] Fix automatic copy even when copy=False in crop --- skimage/util/arraypad.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/util/arraypad.py b/skimage/util/arraypad.py index ee0a3c5d..49f89e3d 100644 --- a/skimage/util/arraypad.py +++ b/skimage/util/arraypad.py @@ -1528,7 +1528,7 @@ def crop(ar, crop_width, copy=False, order='K'): The cropped array. If `copy=False` (default), this is a sliced view of the input array. """ - ar = np.array(ar) + ar = np.array(ar, copy=False) crops = _validate_lengths(ar, crop_width) slices = tuple([slice(a, -b) for a, b in crops]) if copy: From 893c3be289e326cbb70ef0d2cb65847f568fa9c2 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 7 Jan 2015 11:26:21 +1100 Subject: [PATCH 10/36] Import crop into util package namespace --- skimage/util/__init__.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/skimage/util/__init__.py b/skimage/util/__init__.py index 5577e46b..4739bd5b 100644 --- a/skimage/util/__init__.py +++ b/skimage/util/__init__.py @@ -3,7 +3,7 @@ from .dtype import (img_as_float, img_as_int, img_as_uint, img_as_ubyte, from .shape import view_as_blocks, view_as_windows from .noise import random_noise -from .arraypad import pad +from .arraypad import pad, crop from ._regular_grid import regular_grid from .unique import unique_rows @@ -17,6 +17,7 @@ __all__ = ['img_as_float', 'view_as_blocks', 'view_as_windows', 'pad', + 'crop', 'random_noise', 'regular_grid', 'unique_rows'] From 92433b29619783fb1bc064024ca06666b6307dd0 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 7 Jan 2015 11:40:45 +1100 Subject: [PATCH 11/36] Fix and test 0-end bug in crop --- skimage/util/arraypad.py | 2 +- skimage/util/tests/test_arraypad.py | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/skimage/util/arraypad.py b/skimage/util/arraypad.py index 49f89e3d..cce40dc8 100644 --- a/skimage/util/arraypad.py +++ b/skimage/util/arraypad.py @@ -1530,7 +1530,7 @@ def crop(ar, crop_width, copy=False, order='K'): """ ar = np.array(ar, copy=False) crops = _validate_lengths(ar, crop_width) - slices = tuple([slice(a, -b) for a, b in crops]) + slices = [slice(a, ar.shape[i] - b) for i, (a, b) in enumerate(crops)] if copy: cropped = np.copy(ar[slices], order=order) else: diff --git a/skimage/util/tests/test_arraypad.py b/skimage/util/tests/test_arraypad.py index 4d990d0a..d6f993f9 100644 --- a/skimage/util/tests/test_arraypad.py +++ b/skimage/util/tests/test_arraypad.py @@ -1081,5 +1081,11 @@ def test_copy_crop(): assert np.may_share_memory(arr, out1) +def test_zero_crop(): + arr = np.arange(45).reshape(9, 5) + out = crop(arr, 0) + assert out.shape == (9, 5) + + if __name__ == "__main__": np.testing.run_module_suite() From f321fefb8bcf748a5159c63c4fd4ab99880064a2 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 7 Jan 2015 11:41:12 +1100 Subject: [PATCH 12/36] Fix morpho open/close edge bugs by padding image --- skimage/morphology/grey.py | 50 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/skimage/morphology/grey.py b/skimage/morphology/grey.py index c4ff5ffd..ea5f0f07 100644 --- a/skimage/morphology/grey.py +++ b/skimage/morphology/grey.py @@ -1,9 +1,12 @@ """ Grayscale morphological operations """ +import functools + import numpy as np from scipy import ndimage as nd from .misc import default_selem +from ..util import pad, crop __all__ = ['erosion', 'dilation', 'opening', 'closing', 'white_tophat', 'black_tophat'] @@ -79,6 +82,51 @@ def _invert_selem(selem): return inverted +def pad_for_eccentric_selems(func): + """Pad input images for certain morphological operations. + + Parameters + ---------- + func : callable + A morphological function, either opening or closing, that + supports eccentric structuring elements. The inputs must + include at least `image`, `selem`, and `out`. + + Returns + ------- + func_out : callable + A function + + See Also + -------- + ``opening``, ``closing``. + """ + @functools.wraps(func) + def func_out(image, selem, out=None, *args, **kwargs): + pad_widths = [] + padding = False + if out is None: + out = np.empty_like(image) + for axis_len in selem.shape: + if axis_len % 2 == 0: + axis_pad_width = axis_len - 1 + padding = True + else: + axis_pad_width = 0 + pad_widths.append((axis_pad_width,) * 2) + if padding: + image = pad(image, pad_widths, mode='edge') + out_temp = np.empty_like(image) + else: + out_temp = out + out_temp = func(image, selem, out=out_temp, *args, **kwargs) + if padding: + out[:] = crop(out_temp, pad_widths) + else: + out = out_temp + return out + return func_out + @default_selem def erosion(image, selem=None, out=None, shift_x=False, shift_y=False): """Return greyscale morphological erosion of an image. @@ -203,6 +251,7 @@ def dilation(image, selem=None, out=None, shift_x=False, shift_y=False): @default_selem +@pad_for_eccentric_selems def opening(image, selem=None, out=None): """Return greyscale morphological opening of an image. @@ -252,6 +301,7 @@ def opening(image, selem=None, out=None): @default_selem +@pad_for_eccentric_selems def closing(image, selem=None, out=None): """Return greyscale morphological closing of an image. From b1891dc24e8b3d6ab26221dfc521ee9d4d34d57b Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 8 Jan 2015 17:00:06 +1100 Subject: [PATCH 13/36] Make find_boundaries symmetric Fixes #738 --- skimage/segmentation/boundaries.py | 33 +++++++++++--- skimage/segmentation/tests/test_boundaries.py | 43 ++++++++++--------- 2 files changed, 48 insertions(+), 28 deletions(-) diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py index 8e4ef32e..11773b1c 100644 --- a/skimage/segmentation/boundaries.py +++ b/skimage/segmentation/boundaries.py @@ -1,15 +1,34 @@ import numpy as np -from ..morphology import dilation, square +from scipy import ndimage as nd +from ..morphology import dilation, erosion, square from ..util import img_as_float from ..color import gray2rgb from .._shared.utils import deprecated -def find_boundaries(label_img): - """Return bool array where boundaries between labeled regions are True.""" - boundaries = np.zeros(label_img.shape, dtype=np.bool) - boundaries[1:, :] += label_img[1:, :] != label_img[:-1, :] - boundaries[:, 1:] += label_img[:, 1:] != label_img[:, :-1] +def find_boundaries(label_img, connectivity=1): + """Return bool array where boundaries between labeled regions are True. + + Parameters + ---------- + label_img : array of int + An array in which different regions are labeled with different + integers. + connectivity: int in {1, ..., `label_img.ndim`}, optional + A pixel is considered a boundary pixel if any of its neighbors + has a different label. `connectivity` controls which pixels are + considered neighbors. A connectivity of 1 (default) means + pixels sharing an edge (in 2D) or a face (in 3D) will be + considered neighbors. A connectivity of `label_img.ndim` means + pixels sharing a corner will be considered neighbors. + + Returns + ------- + boundaries : array of bool, same shape as `label_img` + A bool image where `True` represents a boundary pixel. + """ + selem = nd.generate_binary_structure(label_img.ndim, connectivity) + boundaries = dilation(label_img, selem) != erosion(label_img, selem) return boundaries @@ -35,7 +54,7 @@ def mark_boundaries(image, label_img, color=(1, 1, 0), boundaries = find_boundaries(label_img) if outline_color is not None: - outer_boundaries = dilation(boundaries.astype(np.uint8), square(2)) + outer_boundaries = dilation(boundaries.astype(np.uint8), square(3)) image[outer_boundaries != 0, :] = np.array(outline_color) image[boundaries, :] = np.array(color) return image diff --git a/skimage/segmentation/tests/test_boundaries.py b/skimage/segmentation/tests/test_boundaries.py index 2fff52f8..7d5e5d9c 100644 --- a/skimage/segmentation/tests/test_boundaries.py +++ b/skimage/segmentation/tests/test_boundaries.py @@ -8,12 +8,12 @@ def test_find_boundaries(): image[2:7, 2:7] = 1 ref = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [0, 0, 1, 1, 1, 1, 1, 1, 0, 0], - [0, 0, 1, 0, 0, 0, 0, 1, 0, 0], - [0, 0, 1, 0, 0, 0, 0, 1, 0, 0], - [0, 0, 1, 0, 0, 0, 0, 1, 0, 0], - [0, 0, 1, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 0, 0], + [0, 1, 1, 0, 0, 0, 1, 1, 0, 0], + [0, 1, 1, 0, 0, 0, 1, 1, 0, 0], + [0, 1, 1, 0, 0, 0, 1, 1, 0, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) @@ -28,27 +28,28 @@ def test_mark_boundaries(): label_image[2:7, 2:7] = 1 ref = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [0, 0, 1, 1, 1, 1, 1, 1, 0, 0], - [0, 0, 1, 0, 0, 0, 0, 1, 0, 0], - [0, 0, 1, 0, 0, 0, 0, 1, 0, 0], - [0, 0, 1, 0, 0, 0, 0, 1, 0, 0], - [0, 0, 1, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 0, 0], + [0, 1, 1, 0, 0, 0, 1, 1, 0, 0], + [0, 1, 1, 0, 0, 0, 1, 1, 0, 0], + [0, 1, 1, 0, 0, 0, 1, 1, 0, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) + result = mark_boundaries(image, label_image, color=(1, 1, 1)).mean(axis=2) assert_array_equal(result, ref) - ref = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [0, 0, 1, 1, 1, 1, 1, 1, 2, 0], - [0, 0, 1, 2, 2, 2, 2, 1, 2, 0], - [0, 0, 1, 2, 0, 0, 0, 1, 2, 0], - [0, 0, 1, 2, 0, 0, 0, 1, 2, 0], - [0, 0, 1, 2, 0, 0, 0, 1, 2, 0], - [0, 0, 1, 1, 1, 1, 1, 2, 2, 0], - [0, 0, 2, 2, 2, 2, 2, 2, 0, 0], + ref = np.array([[0, 2, 2, 2, 2, 2, 2, 2, 0, 0], + [2, 2, 1, 1, 1, 1, 1, 2, 2, 0], + [2, 1, 1, 1, 1, 1, 1, 1, 2, 0], + [2, 1, 1, 2, 2, 2, 1, 1, 2, 0], + [2, 1, 1, 2, 0, 2, 1, 1, 2, 0], + [2, 1, 1, 2, 2, 2, 1, 1, 2, 0], + [2, 1, 1, 1, 1, 1, 1, 1, 2, 0], + [2, 2, 1, 1, 1, 1, 1, 2, 2, 0], + [0, 2, 2, 2, 2, 2, 2, 2, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) result = mark_boundaries(image, label_image, color=(1, 1, 1), outline_color=(2, 2, 2)).mean(axis=2) From d538abdb97470cc4f95c32d6a22acb4f85257de2 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 8 Jan 2015 18:25:51 +1100 Subject: [PATCH 14/36] Fix order kwarg for some versions of numpy.copy It seems it isn't always a valid kwarg. This bypasses that problem by calling `np.array` directly (which is just what `np.copy` does under the hood anyway). --- skimage/util/arraypad.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/util/arraypad.py b/skimage/util/arraypad.py index cce40dc8..3d9d6533 100644 --- a/skimage/util/arraypad.py +++ b/skimage/util/arraypad.py @@ -1532,7 +1532,7 @@ def crop(ar, crop_width, copy=False, order='K'): crops = _validate_lengths(ar, crop_width) slices = [slice(a, ar.shape[i] - b) for i, (a, b) in enumerate(crops)] if copy: - cropped = np.copy(ar[slices], order=order) + cropped = np.array(ar[slices], order=order, copy=True) else: cropped = ar[slices] return cropped From d8555142dab316c082ee86d7c297e7a262a1a7b4 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 8 Jan 2015 22:04:19 +1100 Subject: [PATCH 15/36] Update binary morphology tests, which no longer warn --- skimage/morphology/tests/test_binary.py | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/skimage/morphology/tests/test_binary.py b/skimage/morphology/tests/test_binary.py index d52f92bb..6d8097b4 100644 --- a/skimage/morphology/tests/test_binary.py +++ b/skimage/morphology/tests/test_binary.py @@ -4,7 +4,6 @@ from numpy import testing from skimage import data, color from skimage.util import img_as_bool from skimage.morphology import binary, grey, selem -from skimage._shared._warnings import expected_warnings from scipy import ndimage @@ -15,50 +14,44 @@ bw_img = img > 100 def test_non_square_image(): strel = selem.square(3) binary_res = binary.binary_erosion(bw_img[:100, :200], strel) - with expected_warnings(['precision loss']): - grey_res = img_as_bool(grey.erosion(bw_img[:100, :200], strel)) + grey_res = img_as_bool(grey.erosion(bw_img[:100, :200], strel)) testing.assert_array_equal(binary_res, grey_res) def test_binary_erosion(): strel = selem.square(3) binary_res = binary.binary_erosion(bw_img, strel) - with expected_warnings(['precision loss']): - grey_res = img_as_bool(grey.erosion(bw_img, strel)) + grey_res = img_as_bool(grey.erosion(bw_img, strel)) testing.assert_array_equal(binary_res, grey_res) def test_binary_dilation(): strel = selem.square(3) binary_res = binary.binary_dilation(bw_img, strel) - with expected_warnings(['precision loss']): - grey_res = img_as_bool(grey.dilation(bw_img, strel)) + grey_res = img_as_bool(grey.dilation(bw_img, strel)) testing.assert_array_equal(binary_res, grey_res) def test_binary_closing(): strel = selem.square(3) binary_res = binary.binary_closing(bw_img, strel) - with expected_warnings(['precision loss']): - grey_res = img_as_bool(grey.closing(bw_img, strel)) + grey_res = img_as_bool(grey.closing(bw_img, strel)) testing.assert_array_equal(binary_res, grey_res) def test_binary_opening(): strel = selem.square(3) binary_res = binary.binary_opening(bw_img, strel) - with expected_warnings(['precision loss']): - grey_res = img_as_bool(grey.opening(bw_img, strel)) + grey_res = img_as_bool(grey.opening(bw_img, strel)) testing.assert_array_equal(binary_res, grey_res) def test_selem_overflow(): strel = np.ones((17, 17), dtype=np.uint8) - img = np.zeros((20, 20)) - img[2:19, 2:19] = 1 + img = np.zeros((20, 20), dtype=bool) + img[2:19, 2:19] = True binary_res = binary.binary_erosion(img, strel) - with expected_warnings(['precision loss']): - grey_res = img_as_bool(grey.erosion(img, strel)) + grey_res = img_as_bool(grey.erosion(img, strel)) testing.assert_array_equal(binary_res, grey_res) From ff1631c717cc908248faeb0a32bb6c8ba1527cc7 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 8 Jan 2015 22:05:31 +1100 Subject: [PATCH 16/36] Update find_boundaries with new modes --- skimage/segmentation/boundaries.py | 111 +++++++++++++++++++++++++++-- 1 file changed, 106 insertions(+), 5 deletions(-) diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py index 11773b1c..37fd77fd 100644 --- a/skimage/segmentation/boundaries.py +++ b/skimage/segmentation/boundaries.py @@ -6,7 +6,7 @@ from ..color import gray2rgb from .._shared.utils import deprecated -def find_boundaries(label_img, connectivity=1): +def find_boundaries(label_img, connectivity=1, mode='thick', background=0): """Return bool array where boundaries between labeled regions are True. Parameters @@ -21,15 +21,116 @@ def find_boundaries(label_img, connectivity=1): pixels sharing an edge (in 2D) or a face (in 3D) will be considered neighbors. A connectivity of `label_img.ndim` means pixels sharing a corner will be considered neighbors. + mode: string in {'thick', 'inner', 'outer', 'subpixel'} + How to mark the boundaries: + + - thick: any pixel not completely surrounded by pixels of the + same label (defined by `connectivity`) is marked as a boundary. + This results in boundaries that are 2 pixels thick. + - inner: outline the pixels *just inside* of objects, leaving + background pixels untouched. + - outer: outline pixels in the background around object + boundaries. When two objects touch, their boundary is also + marked. + - subpixel: return a doubled image, with pixels *between* the + original pixels marked as boundary where appropriate. + background: int, optional + For modes 'inner' and 'outer', a definition of a background + label is required. See `mode` for descriptions of these two. Returns ------- boundaries : array of bool, same shape as `label_img` - A bool image where `True` represents a boundary pixel. + A bool image where ``True`` represents a boundary pixel. For + `mode` equal to 'subpixel', ``boundaries.shape[i]`` is equal + to ``2 * label_img.shape[i] - 1`` for all ``i`` (a pixel is + inserted in between all other pairs of pixels). + + Examples + -------- + >>> labels = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ... [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ... [0, 0, 0, 0, 0, 5, 5, 5, 0, 0], + ... [0, 0, 1, 1, 1, 5, 5, 5, 0, 0], + ... [0, 0, 1, 1, 1, 5, 5, 5, 0, 0], + ... [0, 0, 1, 1, 1, 5, 5, 5, 0, 0], + ... [0, 0, 0, 0, 0, 5, 5, 5, 0, 0], + ... [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + ... [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=np.uint8) + >>> find_boundaries(labels).astype(np.uint8) # display 1/0, not True/False + array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 1, 1, 0, 0], + [0, 0, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 0, 1, 1, 0], + [0, 1, 1, 0, 1, 1, 0, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 0, 1, 1, 0], + [0, 0, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 0, 0, 0, 0, 1, 1, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) + >>> find_boundaries(labels, mode='inner').astype(np.uint8) + array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 1, 1, 0, 0], + [0, 0, 1, 1, 1, 1, 0, 1, 0, 0], + [0, 0, 1, 0, 1, 1, 0, 1, 0, 0], + [0, 0, 1, 1, 1, 1, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 1, 1, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) + >>> find_boundaries(labels, mode='outer').astype(np.uint8) + array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 1, 1, 0, 0], + [0, 0, 1, 1, 1, 1, 0, 0, 1, 0], + [0, 1, 0, 0, 1, 1, 0, 0, 1, 0], + [0, 1, 0, 0, 1, 1, 0, 0, 1, 0], + [0, 1, 0, 0, 1, 1, 0, 0, 1, 0], + [0, 0, 1, 1, 1, 1, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 1, 1, 1, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) + >>> find_boundaries(labels, mode='subpixel').astype(np.uint8) + array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) """ - selem = nd.generate_binary_structure(label_img.ndim, connectivity) - boundaries = dilation(label_img, selem) != erosion(label_img, selem) - return boundaries + ndim = label_img.ndim + selem = nd.generate_binary_structure(ndim, connectivity) + if mode != 'subpixel': + boundaries = dilation(label_img, selem) != erosion(label_img, selem) + if mode == 'inner': + foreground_image = (label_img != background) + boundaries &= foreground_image + elif mode == 'outer': + background_image = (label_img == background) + selem = nd.generate_binary_structure(ndim, ndim) + no_adjacent_background = ~dilation(background_image, selem) + boundaries &= (background_image | no_adjacent_background) + return boundaries + else: + label_img_expanded = np.zeros([(2 * s - 1) for s in label_img.shape], + label_img.dtype) + pixels = [slice(None, None, 2)] * ndim + selem = nd.generate_binary_structure(ndim, ndim) + label_img_expanded[pixels] = label_img + max_label = np.iinfo(label_img.dtype).max + label_img_edge_inverted = np.array(label_img_expanded, copy=True) + label_img_edge_inverted[label_img_expanded == 0] = max_label + boundaries = (dilation(label_img_expanded, selem) != + erosion(label_img_edge_inverted, selem)) + return boundaries def mark_boundaries(image, label_img, color=(1, 1, 0), From 32a960db0d3d5e0f93d1cee27961043bd479f50c Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 8 Jan 2015 22:06:31 +1100 Subject: [PATCH 17/36] Fix 'outer' mode of find_boundaries --- skimage/segmentation/boundaries.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py index 37fd77fd..89127401 100644 --- a/skimage/segmentation/boundaries.py +++ b/skimage/segmentation/boundaries.py @@ -108,6 +108,7 @@ def find_boundaries(label_img, connectivity=1, mode='thick', background=0): """ ndim = label_img.ndim selem = nd.generate_binary_structure(ndim, connectivity) + max_label = np.iinfo(label_img.dtype).max if mode != 'subpixel': boundaries = dilation(label_img, selem) != erosion(label_img, selem) if mode == 'inner': @@ -116,8 +117,12 @@ def find_boundaries(label_img, connectivity=1, mode='thick', background=0): elif mode == 'outer': background_image = (label_img == background) selem = nd.generate_binary_structure(ndim, ndim) - no_adjacent_background = ~dilation(background_image, selem) - boundaries &= (background_image | no_adjacent_background) + inverted_background = np.array(label_img, copy=True) + inverted_background[background_image] = max_label + adjacent_objects = ((dilation(label_img, selem) != + erosion(inverted_background, selem)) & + ~background_image) + boundaries &= (background_image | adjacent_objects) return boundaries else: label_img_expanded = np.zeros([(2 * s - 1) for s in label_img.shape], From 9c78fce4cd253162a06613384b49c3ad0428af9c Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Sun, 11 Jan 2015 11:51:29 +1100 Subject: [PATCH 18/36] Remove unused deprecation warning --- skimage/segmentation/boundaries.py | 1 - 1 file changed, 1 deletion(-) diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py index 89127401..255f9816 100644 --- a/skimage/segmentation/boundaries.py +++ b/skimage/segmentation/boundaries.py @@ -3,7 +3,6 @@ from scipy import ndimage as nd from ..morphology import dilation, erosion, square from ..util import img_as_float from ..color import gray2rgb -from .._shared.utils import deprecated def find_boundaries(label_img, connectivity=1, mode='thick', background=0): From fb893c277d572ec19d9b8f44d2214c5b1bbe9fbd Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Sun, 11 Jan 2015 18:53:42 +1100 Subject: [PATCH 19/36] Add correct subpixel boundary estimation. --- skimage/segmentation/boundaries.py | 52 +++++++++++++++++++++++------- 1 file changed, 41 insertions(+), 11 deletions(-) diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py index 255f9816..306e213d 100644 --- a/skimage/segmentation/boundaries.py +++ b/skimage/segmentation/boundaries.py @@ -1,10 +1,49 @@ import numpy as np from scipy import ndimage as nd from ..morphology import dilation, erosion, square -from ..util import img_as_float +from ..util import img_as_float, view_as_windows, pad from ..color import gray2rgb +def _find_boundaries_subpixel(label_img): + """See ``find_boundaries(..., mode='subpixel')``. + + Notes + ----- + This function puts in an empty row and column between each *actual* + row and column of the image, for a corresponding shape of $2s - 1$ + for every image dimension of size $s$. These "interstitial" rows + and columns are filled as ``True`` if they separate two labels in + `label_img`, ``False`` otherwise. + + I used ``view_as_windows`` to get the neighborhood of each pixel. + Then I check whether there are two labels or more in that + neighborhood. + """ + ndim = label_img.ndim + max_label = np.iinfo(label_img.dtype).max + + label_img_expanded = np.zeros([(2 * s - 1) for s in label_img.shape], + label_img.dtype) + pixels = [slice(None, None, 2)] * ndim + label_img_expanded[pixels] = label_img + + edges = np.ones(label_img_expanded.shape, dtype=bool) + edges[pixels] = False + label_img_expanded[edges] = max_label + windows = view_as_windows(pad(label_img_expanded, 1, + mode='constant', constant_values=0), + (3,) * ndim) + + boundaries = np.zeros_like(edges) + for index in np.ndindex(label_img_expanded.shape): + if edges[index]: + values = np.unique(windows[index].ravel()) + if len(values) > 2: # single value and max_label + boundaries[index] = True + return boundaries + + def find_boundaries(label_img, connectivity=1, mode='thick', background=0): """Return bool array where boundaries between labeled regions are True. @@ -124,16 +163,7 @@ def find_boundaries(label_img, connectivity=1, mode='thick', background=0): boundaries &= (background_image | adjacent_objects) return boundaries else: - label_img_expanded = np.zeros([(2 * s - 1) for s in label_img.shape], - label_img.dtype) - pixels = [slice(None, None, 2)] * ndim - selem = nd.generate_binary_structure(ndim, ndim) - label_img_expanded[pixels] = label_img - max_label = np.iinfo(label_img.dtype).max - label_img_edge_inverted = np.array(label_img_expanded, copy=True) - label_img_edge_inverted[label_img_expanded == 0] = max_label - boundaries = (dilation(label_img_expanded, selem) != - erosion(label_img_edge_inverted, selem)) + boundaries = _find_boundaries_subpixel(label_img) return boundaries From a55e3346f7bfc7e836126eeb34bbaa6efb8b5d43 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Sun, 11 Jan 2015 19:22:32 +1100 Subject: [PATCH 20/36] Set up correct type for find_boundaries test --- skimage/segmentation/tests/test_boundaries.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/segmentation/tests/test_boundaries.py b/skimage/segmentation/tests/test_boundaries.py index 7d5e5d9c..4019e112 100644 --- a/skimage/segmentation/tests/test_boundaries.py +++ b/skimage/segmentation/tests/test_boundaries.py @@ -4,7 +4,7 @@ from skimage.segmentation import find_boundaries, mark_boundaries def test_find_boundaries(): - image = np.zeros((10, 10)) + image = np.zeros((10, 10), dtype=np.uint8) image[2:7, 2:7] = 1 ref = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], @@ -24,7 +24,7 @@ def test_find_boundaries(): def test_mark_boundaries(): image = np.zeros((10, 10)) - label_image = np.zeros((10, 10)) + label_image = np.zeros((10, 10), dtype=np.uint8) label_image[2:7, 2:7] = 1 ref = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], From a406f8dd2a5065f39af300e47b3bf3e60ef1565b Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Sun, 11 Jan 2015 19:23:17 +1100 Subject: [PATCH 21/36] Make subpixel find_boundary doctest more readable --- skimage/segmentation/boundaries.py | 37 +++++++++++++++--------------- 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py index 306e213d..f4cca6d5 100644 --- a/skimage/segmentation/boundaries.py +++ b/skimage/segmentation/boundaries.py @@ -125,34 +125,33 @@ def find_boundaries(label_img, connectivity=1, mode='thick', background=0): [0, 0, 1, 1, 1, 1, 0, 0, 1, 0], [0, 0, 0, 0, 0, 1, 1, 1, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) - >>> find_boundaries(labels, mode='subpixel').astype(np.uint8) - array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], - [0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], - [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], - [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], - [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], - [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], - [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], - [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], - [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], - [0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) + >>> labels_small = labels[::2, ::3] + >>> labels_small + array([[0, 0, 0, 0], + [0, 0, 5, 0], + [0, 1, 5, 0], + [0, 0, 5, 0], + [0, 0, 0, 0]], dtype=uint8) + >>> find_boundaries(labels_small, mode='subpixel').astype(np.uint8) + array([[0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 1, 1, 1, 0], + [0, 0, 0, 1, 0, 1, 0], + [0, 1, 1, 1, 0, 1, 0], + [0, 1, 0, 1, 0, 1, 0], + [0, 1, 1, 1, 0, 1, 0], + [0, 0, 0, 1, 0, 1, 0], + [0, 0, 0, 1, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 0]], dtype=uint8) """ ndim = label_img.ndim selem = nd.generate_binary_structure(ndim, connectivity) - max_label = np.iinfo(label_img.dtype).max if mode != 'subpixel': boundaries = dilation(label_img, selem) != erosion(label_img, selem) if mode == 'inner': foreground_image = (label_img != background) boundaries &= foreground_image elif mode == 'outer': + max_label = np.iinfo(label_img.dtype).max background_image = (label_img == background) selem = nd.generate_binary_structure(ndim, ndim) inverted_background = np.array(label_img, copy=True) From 28591cc3ee070ede16fd17a082a102f8f3f305a6 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 12 Jan 2015 12:08:32 +1100 Subject: [PATCH 22/36] Allow different modes in mark_boundaries --- skimage/segmentation/boundaries.py | 38 +++++++++++++++++++++--------- 1 file changed, 27 insertions(+), 11 deletions(-) diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py index f4cca6d5..583d593b 100644 --- a/skimage/segmentation/boundaries.py +++ b/skimage/segmentation/boundaries.py @@ -1,3 +1,5 @@ +from __future__ import division + import numpy as np from scipy import ndimage as nd from ..morphology import dilation, erosion, square @@ -95,7 +97,7 @@ def find_boundaries(label_img, connectivity=1, mode='thick', background=0): ... [0, 0, 0, 0, 0, 5, 5, 5, 0, 0], ... [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], ... [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=np.uint8) - >>> find_boundaries(labels).astype(np.uint8) # display 1/0, not True/False + >>> find_boundaries(labels, mode='thick').astype(np.uint8) array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 1, 1, 0], @@ -167,28 +169,42 @@ def find_boundaries(label_img, connectivity=1, mode='thick', background=0): def mark_boundaries(image, label_img, color=(1, 1, 0), - outline_color=(0, 0, 0)): + outline_color=(0, 0, 0), mode='outer'): """Return image with boundaries between labeled regions highlighted. Parameters ---------- image : (M, N[, 3]) array Grayscale or RGB image. - label_img : (M, N) array + label_img : (M, N) array of int Label array where regions are marked by different integer values. - color : length-3 sequence + color : length-3 sequence, optional RGB color of boundaries in the output image. - outline_color : length-3 sequence + outline_color : length-3 sequence, optional RGB color surrounding boundaries in the output image. If None, no outline is drawn. + mode : string in {'thick', 'inner', 'outer', 'subpixel'}, optional + The mode for finding boundaries. + + Returns + ------- + marked : (M, N, 3) array of float + An image in which the boundaries between labels are + superimposed on the original image. + + See Also + -------- + ``find_boundaries``. """ if image.ndim == 2: image = gray2rgb(image) - image = img_as_float(image, force_copy=True) - - boundaries = find_boundaries(label_img) + marked = img_as_float(image, force_copy=True) + if mode == 'subpixel': + marked = nd.zoom(marked, [2 - 1/s for s in marked.shape[:-1]] + [1], + mode='reflect') + boundaries = find_boundaries(label_img, mode=mode) if outline_color is not None: outer_boundaries = dilation(boundaries.astype(np.uint8), square(3)) - image[outer_boundaries != 0, :] = np.array(outline_color) - image[boundaries, :] = np.array(color) - return image + marked[outer_boundaries != 0, :] = np.array(outline_color) + marked[boundaries, :] = np.array(color) + return marked From 257502c099d77e0d3bac2ed24dde5ed093a0345d Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 12 Jan 2015 12:09:22 +1100 Subject: [PATCH 23/36] Prevent unnecessary array copy for RGB images --- skimage/segmentation/boundaries.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py index 583d593b..67763bbb 100644 --- a/skimage/segmentation/boundaries.py +++ b/skimage/segmentation/boundaries.py @@ -196,9 +196,9 @@ def mark_boundaries(image, label_img, color=(1, 1, 0), -------- ``find_boundaries``. """ - if image.ndim == 2: - image = gray2rgb(image) marked = img_as_float(image, force_copy=True) + if marked.ndim == 2: + marked = gray2rgb(marked) if mode == 'subpixel': marked = nd.zoom(marked, [2 - 1/s for s in marked.shape[:-1]] + [1], mode='reflect') From cb2002c36777be701c37d1958d5218ce289403ae Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 12 Jan 2015 12:15:33 +1100 Subject: [PATCH 24/36] Allow custom background_label in mark_boundaries --- skimage/segmentation/boundaries.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py index 67763bbb..49764432 100644 --- a/skimage/segmentation/boundaries.py +++ b/skimage/segmentation/boundaries.py @@ -169,7 +169,7 @@ def find_boundaries(label_img, connectivity=1, mode='thick', background=0): def mark_boundaries(image, label_img, color=(1, 1, 0), - outline_color=(0, 0, 0), mode='outer'): + outline_color=(0, 0, 0), mode='outer', background_label=0): """Return image with boundaries between labeled regions highlighted. Parameters @@ -185,6 +185,9 @@ def mark_boundaries(image, label_img, color=(1, 1, 0), outline is drawn. mode : string in {'thick', 'inner', 'outer', 'subpixel'}, optional The mode for finding boundaries. + background_label : int, optional + Which label to consider background (this is only useful for + modes ``inner`` and ``outer``). Returns ------- @@ -202,7 +205,8 @@ def mark_boundaries(image, label_img, color=(1, 1, 0), if mode == 'subpixel': marked = nd.zoom(marked, [2 - 1/s for s in marked.shape[:-1]] + [1], mode='reflect') - boundaries = find_boundaries(label_img, mode=mode) + boundaries = find_boundaries(label_img, mode=mode, + background=background_label) if outline_color is not None: outer_boundaries = dilation(boundaries.astype(np.uint8), square(3)) marked[outer_boundaries != 0, :] = np.array(outline_color) From 242dbf9b2bc02cddda942c1f0a0b03151eb232f2 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 12 Jan 2015 12:24:17 +1100 Subject: [PATCH 25/36] Don't outline boundaries by default mark_boundaries used to draw outlines by default. This was necessary partly because of the asymmetric way boundaries were being drawn. Now, boundaries are either thick or drawn around background or drawn with subpixel precision, so I suggest that outlines are less needed. --- skimage/segmentation/boundaries.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py index 49764432..c637dd5f 100644 --- a/skimage/segmentation/boundaries.py +++ b/skimage/segmentation/boundaries.py @@ -169,7 +169,7 @@ def find_boundaries(label_img, connectivity=1, mode='thick', background=0): def mark_boundaries(image, label_img, color=(1, 1, 0), - outline_color=(0, 0, 0), mode='outer', background_label=0): + outline_color=None, mode='outer', background_label=0): """Return image with boundaries between labeled regions highlighted. Parameters From 307d4b6137962910c7fffdad32e10ee37b3324e2 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 12 Jan 2015 12:26:41 +1100 Subject: [PATCH 26/36] Simplify code for coloring boundaries --- skimage/segmentation/boundaries.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py index c637dd5f..a2c5f3fb 100644 --- a/skimage/segmentation/boundaries.py +++ b/skimage/segmentation/boundaries.py @@ -208,7 +208,7 @@ def mark_boundaries(image, label_img, color=(1, 1, 0), boundaries = find_boundaries(label_img, mode=mode, background=background_label) if outline_color is not None: - outer_boundaries = dilation(boundaries.astype(np.uint8), square(3)) - marked[outer_boundaries != 0, :] = np.array(outline_color) - marked[boundaries, :] = np.array(color) + outlines = dilation(boundaries, square(3)) + marked[outlines] = outline_color + marked[boundaries] = color return marked From 6586d1eb3fd01367cad5b73b0430a2a9bcfe23c7 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 12 Jan 2015 15:31:56 +1100 Subject: [PATCH 27/36] Test subpixel mark_boundaries --- skimage/segmentation/tests/test_boundaries.py | 33 +++++++++++++++++-- 1 file changed, 30 insertions(+), 3 deletions(-) diff --git a/skimage/segmentation/tests/test_boundaries.py b/skimage/segmentation/tests/test_boundaries.py index 4019e112..7b837f4e 100644 --- a/skimage/segmentation/tests/test_boundaries.py +++ b/skimage/segmentation/tests/test_boundaries.py @@ -1,8 +1,11 @@ import numpy as np -from numpy.testing import assert_array_equal +from numpy.testing import assert_array_equal, assert_allclose from skimage.segmentation import find_boundaries, mark_boundaries +white = (1, 1, 1) + + def test_find_boundaries(): image = np.zeros((10, 10), dtype=np.uint8) image[2:7, 2:7] = 1 @@ -38,7 +41,7 @@ def test_mark_boundaries(): [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) - result = mark_boundaries(image, label_image, color=(1, 1, 1)).mean(axis=2) + result = mark_boundaries(image, label_image, color=white).mean(axis=2) assert_array_equal(result, ref) ref = np.array([[0, 2, 2, 2, 2, 2, 2, 2, 0, 0], @@ -51,10 +54,34 @@ def test_mark_boundaries(): [2, 2, 1, 1, 1, 1, 1, 2, 2, 0], [0, 2, 2, 2, 2, 2, 2, 2, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) - result = mark_boundaries(image, label_image, color=(1, 1, 1), + result = mark_boundaries(image, label_image, color=white, outline_color=(2, 2, 2)).mean(axis=2) assert_array_equal(result, ref) +def test_mark_boundaries_subpixel(): + labels = np.array([[0, 0, 0, 0], + [0, 0, 5, 0], + [0, 1, 5, 0], + [0, 0, 5, 0], + [0, 0, 0, 0]], dtype=np.uint8) + np.random.seed(0) + image = np.round(np.random.rand(*labels.shape), 2) + marked = np.mark_boundaries(image, labels, color=white, mode='subpixel') + marked_proj = np.round(np.mean(marked, axis=-1), 2) + + ref_result = np.array( + [[ 0.55, 0.63, 0.72, 0.69, 0.6 , 0.55, 0.54], + [ 0.45, 0.58, 0.72, 1. , 1. , 1. , 0.69], + [ 0.42, 0.54, 0.65, 1. , 0.44, 1. , 0.89], + [ 0.69, 1. , 1. , 1. , 0.69, 1. , 0.83], + [ 0.96, 1. , 0.38, 1. , 0.79, 1. , 0.53], + [ 0.89, 1. , 1. , 1. , 0.38, 1. , 0.16], + [ 0.57, 0.78, 0.93, 1. , 0.07, 1. , 0.09], + [ 0.2 , 0.52, 0.92, 1. , 1. , 1. , 0.54], + [ 0.02, 0.35, 0.83, 0.9 , 0.78, 0.81, 0.87]]) + assert_allclose(marked_proj, ref_result, atol=0.01) + + if __name__ == "__main__": np.testing.run_module_suite() From 9b9730b058c8f1b11cfcf2fe1438b2c672874e28 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 12 Jan 2015 15:59:38 +1100 Subject: [PATCH 28/36] Specify mode for mark_boundary tests --- skimage/segmentation/tests/test_boundaries.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/skimage/segmentation/tests/test_boundaries.py b/skimage/segmentation/tests/test_boundaries.py index 7b837f4e..0992fa74 100644 --- a/skimage/segmentation/tests/test_boundaries.py +++ b/skimage/segmentation/tests/test_boundaries.py @@ -41,7 +41,8 @@ def test_mark_boundaries(): [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) - result = mark_boundaries(image, label_image, color=white).mean(axis=2) + marked = mark_boundaries(image, label_image, color=white, mode='thick') + result = np.mean(marked, axis=-1) assert_array_equal(result, ref) ref = np.array([[0, 2, 2, 2, 2, 2, 2, 2, 0, 0], @@ -54,8 +55,9 @@ def test_mark_boundaries(): [2, 2, 1, 1, 1, 1, 1, 2, 2, 0], [0, 2, 2, 2, 2, 2, 2, 2, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]) - result = mark_boundaries(image, label_image, color=white, - outline_color=(2, 2, 2)).mean(axis=2) + marked = mark_boundaries(image, label_image, color=white, + outline_color=(2, 2, 2), mode='thick') + result = np.mean(marked, axis=-1) assert_array_equal(result, ref) From ac5c0c30a985432e34cbb9d69c43221bb88209c3 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Mon, 12 Jan 2015 16:01:01 +1100 Subject: [PATCH 29/36] Fix incorrect call to mark_boundaries in test --- skimage/segmentation/tests/test_boundaries.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/segmentation/tests/test_boundaries.py b/skimage/segmentation/tests/test_boundaries.py index 0992fa74..e6e401ac 100644 --- a/skimage/segmentation/tests/test_boundaries.py +++ b/skimage/segmentation/tests/test_boundaries.py @@ -69,7 +69,7 @@ def test_mark_boundaries_subpixel(): [0, 0, 0, 0]], dtype=np.uint8) np.random.seed(0) image = np.round(np.random.rand(*labels.shape), 2) - marked = np.mark_boundaries(image, labels, color=white, mode='subpixel') + marked = mark_boundaries(image, labels, color=white, mode='subpixel') marked_proj = np.round(np.mean(marked, axis=-1), 2) ref_result = np.array( From 7935e515968246e2b07d5e0537eaa51dfdc7bfce Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 14 Jan 2015 17:24:43 +1100 Subject: [PATCH 30/36] Clarify and fix docstrings --- skimage/morphology/binary.py | 2 +- skimage/morphology/grey.py | 27 ++++++++++++++------------- skimage/segmentation/boundaries.py | 4 ++++ skimage/util/arraypad.py | 22 +++++++++++----------- 4 files changed, 30 insertions(+), 25 deletions(-) diff --git a/skimage/morphology/binary.py b/skimage/morphology/binary.py index b12b2948..0ff5c8f7 100644 --- a/skimage/morphology/binary.py +++ b/skimage/morphology/binary.py @@ -7,7 +7,7 @@ from .misc import default_selem # The default_selem decorator provides a diamond structuring element as default -# with the appropriate dimension for the input `image`. +# with the same dimension as the input image and size 3 along each axis. @default_selem def binary_erosion(image, selem=None, out=None): """Return fast binary morphological erosion of an image. diff --git a/skimage/morphology/grey.py b/skimage/morphology/grey.py index ea5f0f07..dfe7dc48 100644 --- a/skimage/morphology/grey.py +++ b/skimage/morphology/grey.py @@ -15,8 +15,8 @@ __all__ = ['erosion', 'dilation', 'opening', 'closing', 'white_tophat', def _shift_selem(selem, shift_x, shift_y): """Shift the binary image `selem` in the left and/or up. - This only affects structuring elements with even number of rows or - columns. + This only affects 2D structuring elements with even number of rows + or columns. Parameters ---------- @@ -63,7 +63,7 @@ def _invert_selem(selem): Returns ------- - inverted : array + inverted : array, same shape and type as `selem` The structuring element, in opposite order. Examples @@ -89,13 +89,14 @@ def pad_for_eccentric_selems(func): ---------- func : callable A morphological function, either opening or closing, that - supports eccentric structuring elements. The inputs must + supports eccentric structuring elements. Its parameters must include at least `image`, `selem`, and `out`. Returns ------- func_out : callable - A function + The same function, but correctly padding the input image before + applying the input function. See Also -------- @@ -156,9 +157,9 @@ def erosion(image, selem=None, out=None, shift_x=False, shift_y=False): Notes ----- - For `uint8` (and `uint16` up to a certain bit-depth) data, the lower - algorithm complexity makes the `skimage.filter.rank.minimum` function more - efficient for larger images and structuring elements. + For ``uint8`` (and ``uint16`` up to a certain bit-depth) data, the + lower algorithm complexity makes the `skimage.filter.rank.minimum` + function more efficient for larger images and structuring elements. Examples -------- @@ -211,7 +212,7 @@ def dilation(image, selem=None, out=None, shift_x=False, shift_y=False): Returns ------- - dilated : uint8 array + dilated : uint8 array, same shape and type as `image` The result of the morphological dilation. Notes @@ -273,7 +274,7 @@ def opening(image, selem=None, out=None): Returns ------- - opening : array + opening : array, same shape and type as `image` The result of the morphological opening. Examples @@ -323,7 +324,7 @@ def closing(image, selem=None, out=None): Returns ------- - closing : array + closing : array, same shape and type as `image` The result of the morphological closing. Examples @@ -371,7 +372,7 @@ def white_tophat(image, selem=None, out=None): Returns ------- - out : array + out : array, same shape and type as `image` The result of the morphological white top hat. Examples @@ -425,7 +426,7 @@ def black_tophat(image, selem=None, out=None): Returns ------- - opening : uint8 array + opening : array, same shape and type as `image` The result of the black top filter. Examples diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py index a2c5f3fb..28cc0279 100644 --- a/skimage/segmentation/boundaries.py +++ b/skimage/segmentation/boundaries.py @@ -203,6 +203,10 @@ def mark_boundaries(image, label_img, color=(1, 1, 0), if marked.ndim == 2: marked = gray2rgb(marked) if mode == 'subpixel': + # Here, we want to interpose an extra line of pixels between + # each original line — except for the last axis which holds + # the RGB information. ``nd.zoom`` then performs the (cubic) + # interpolation, filling in the values of the interposed pixels marked = nd.zoom(marked, [2 - 1/s for s in marked.shape[:-1]] + [1], mode='reflect') boundaries = find_boundaries(label_img, mode=mode, diff --git a/skimage/util/arraypad.py b/skimage/util/arraypad.py index 3d9d6533..f8104df5 100644 --- a/skimage/util/arraypad.py +++ b/skimage/util/arraypad.py @@ -1507,25 +1507,25 @@ def crop(ar, crop_width, copy=False, order='K'): Input array. crop_width : {sequence, int} Number of values to remove from the edges of each axis. - ((before_1, after_1), ... (before_N, after_N)) specifies unique - crop widths at the start and end of each axis. - ((before, after),) specifies the same start and end crops for - every axis. - (int,) or int is a shortcut for before = after = int for all - axes. + ``((before_1, after_1),`` ... ``(before_N, after_N))`` specifies + unique crop widths at the start and end of each axis. + ``((before, after),)`` specifies a fixed start and end crop + for every axis. + ``(n,)`` or ``n`` for integer ``n`` is a shortcut for + before = after = ``n`` for all axes. copy : bool, optional - Ensure that the returned array is contiguous. Normally, a crop + Ensure the returned array is a contiguous copy. Normally, a crop operation will return a discontiguous view of the underlying - input array. Passing `copy=True` will result in a contiguous + input array. Passing ``copy=True`` will result in a contiguous copy. order : {'C', 'F', 'A', 'K'}, optional - If `copy==True`, control the memory layout of the copy. See - `np.copy`. + If ``copy==True``, control the memory layout of the copy. See + ``np.copy``. Returns ------- cropped : array - The cropped array. If `copy=False` (default), this is a sliced + The cropped array. If ``copy=False`` (default), this is a sliced view of the input array. """ ar = np.array(ar, copy=False) From 9eed552ee79d67c4e86b65d154d83b85f9ec48fc Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Wed, 14 Jan 2015 17:50:31 +1100 Subject: [PATCH 31/36] Replace em-dash with regular dash in comment --- skimage/segmentation/boundaries.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py index 28cc0279..1975d879 100644 --- a/skimage/segmentation/boundaries.py +++ b/skimage/segmentation/boundaries.py @@ -204,7 +204,7 @@ def mark_boundaries(image, label_img, color=(1, 1, 0), marked = gray2rgb(marked) if mode == 'subpixel': # Here, we want to interpose an extra line of pixels between - # each original line — except for the last axis which holds + # each original line - except for the last axis which holds # the RGB information. ``nd.zoom`` then performs the (cubic) # interpolation, filling in the values of the interposed pixels marked = nd.zoom(marked, [2 - 1/s for s in marked.shape[:-1]] + [1], From 2d0a3eff339f19322f493a418be8818616ac5c35 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 22 Jan 2015 09:26:06 +1100 Subject: [PATCH 32/36] Expand explanation of _invert_selem --- skimage/morphology/grey.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/skimage/morphology/grey.py b/skimage/morphology/grey.py index dfe7dc48..1a59b934 100644 --- a/skimage/morphology/grey.py +++ b/skimage/morphology/grey.py @@ -241,9 +241,11 @@ def dilation(image, selem=None, out=None, shift_x=False, shift_y=False): """ selem = np.array(selem) selem = _shift_selem(selem, shift_x, shift_y) - # invert the structuring element to patch the same thing happening inside - # ndimage.grey_dilation - # https://github.com/scipy/scipy/blob/ec20ababa400e39ac3ffc9148c01ef86d5349332/scipy/ndimage/morphology.py#L1285 + # Inside ndimage.grey_dilation, the structuring element is inverted, + # eg. `selem = selem[::-1, ::-1]` for 2D [1]_, for reasons unknown to + # this author (@jni). To "patch" this behaviour, we invert our own + # selem before passing it to `nd.grey_dilation`. + # [1] https://github.com/scipy/scipy/blob/ec20ababa400e39ac3ffc9148c01ef86d5349332/scipy/ndimage/morphology.py#L1285 selem = _invert_selem(selem) if out is None: out = np.empty_like(image) From 658a83d46c08e676de5904ef1f26172789d5fe3d Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 22 Jan 2015 11:18:38 +1100 Subject: [PATCH 33/36] Don't use backticks to refer to functions See the "See Also" section in the NumPy Docstring Conventions document: https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt --- skimage/morphology/grey.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/morphology/grey.py b/skimage/morphology/grey.py index 1a59b934..8ab162cf 100644 --- a/skimage/morphology/grey.py +++ b/skimage/morphology/grey.py @@ -100,7 +100,7 @@ def pad_for_eccentric_selems(func): See Also -------- - ``opening``, ``closing``. + opening, closing. """ @functools.wraps(func) def func_out(image, selem, out=None, *args, **kwargs): From 1439341461de15a879aa280c6a73ffebf7b5157d Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 22 Jan 2015 11:27:23 +1100 Subject: [PATCH 34/36] Remove unused cmorph extension The C morphology module has been superseded by wrapping scipy.ndimage. --- bento.info | 3 --- skimage/filters/rank/README.rst | 2 +- skimage/filters/rank/tests/test_rank.py | 10 +++++----- skimage/morphology/setup.py | 3 --- 4 files changed, 6 insertions(+), 12 deletions(-) diff --git a/bento.info b/bento.info index 6a8a0000..a36a4ebb 100644 --- a/bento.info +++ b/bento.info @@ -79,9 +79,6 @@ Library: Extension: skimage.graph._spath Sources: skimage/graph/_spath.pyx - Extension: skimage.morphology.cmorph - Sources: - skimage/morphology/cmorph.pyx Extension: skimage.graph.heap Sources: skimage/graph/heap.pyx diff --git a/skimage/filters/rank/README.rst b/skimage/filters/rank/README.rst index e5c5a9ad..2d6bfb25 100644 --- a/skimage/filters/rank/README.rst +++ b/skimage/filters/rank/README.rst @@ -16,7 +16,7 @@ followed by the moving window is given hereunder /--------------------------/ \-------------------------- ... -We compare cmorph.dilate to this histogram based method to show how +We compare grey.dilate to this histogram based method to show how computational costs increase with respect to image size or structuring element size. This implementation gives better results for large structuring elements. diff --git a/skimage/filters/rank/tests/test_rank.py b/skimage/filters/rank/tests/test_rank.py index 64c3ca97..0fda186c 100644 --- a/skimage/filters/rank/tests/test_rank.py +++ b/skimage/filters/rank/tests/test_rank.py @@ -5,7 +5,7 @@ from numpy.testing import run_module_suite, assert_equal, assert_raises import skimage from skimage import img_as_ubyte, img_as_float from skimage import data, util, morphology -from skimage.morphology import cmorph, disk +from skimage.morphology import grey, disk from skimage.filters import rank from skimage._shared._warnings import expected_warnings @@ -118,7 +118,7 @@ def test_random_sizes(): assert_equal(image16.shape, out16.shape) -def test_compare_with_cmorph_dilate(): +def test_compare_with_grey_dilation(): # compare the result of maximum filter with dilate image = (np.random.rand(100, 100) * 256).astype(np.uint8) @@ -128,11 +128,11 @@ def test_compare_with_cmorph_dilate(): for r in range(1, 20, 1): elem = np.ones((r, r), dtype=np.uint8) rank.maximum(image=image, selem=elem, out=out, mask=mask) - cm = cmorph._dilate(image=image, selem=elem) + cm = grey.dilation(image=image, selem=elem) assert_equal(out, cm) -def test_compare_with_cmorph_erode(): +def test_compare_with_grey_erosion(): # compare the result of maximum filter with erode image = (np.random.rand(100, 100) * 256).astype(np.uint8) @@ -142,7 +142,7 @@ def test_compare_with_cmorph_erode(): for r in range(1, 20, 1): elem = np.ones((r, r), dtype=np.uint8) rank.minimum(image=image, selem=elem, out=out, mask=mask) - cm = cmorph._erode(image=image, selem=elem) + cm = grey.erosion(image=image, selem=elem) assert_equal(out, cm) diff --git a/skimage/morphology/setup.py b/skimage/morphology/setup.py index de1a8794..dbbcad8b 100644 --- a/skimage/morphology/setup.py +++ b/skimage/morphology/setup.py @@ -12,14 +12,11 @@ def configuration(parent_package='', top_path=None): config = Configuration('morphology', parent_package, top_path) config.add_data_dir('tests') - cython(['cmorph.pyx'], working_path=base_path) cython(['_watershed.pyx'], working_path=base_path) cython(['_skeletonize_cy.pyx'], working_path=base_path) cython(['_convex_hull.pyx'], working_path=base_path) cython(['_greyreconstruct.pyx'], working_path=base_path) - config.add_extension('cmorph', sources=['cmorph.c'], - include_dirs=[get_numpy_include_dirs()]) config.add_extension('_watershed', sources=['_watershed.c'], include_dirs=[get_numpy_include_dirs()]) config.add_extension('_skeletonize_cy', sources=['_skeletonize_cy.c'], From ea989e2456d8e9aee7521a410f475ccb1f53ebdb Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 22 Jan 2015 11:28:13 +1100 Subject: [PATCH 35/36] Fix a few typos --- skimage/filters/rank/README.rst | 2 +- skimage/filters/rank/tests/test_rank.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/skimage/filters/rank/README.rst b/skimage/filters/rank/README.rst index 2d6bfb25..6e063ec2 100644 --- a/skimage/filters/rank/README.rst +++ b/skimage/filters/rank/README.rst @@ -26,7 +26,7 @@ update the local histogram. The histogram size is 8-bit (256 bins) for 8-bit images and 2 to 16-bit for 16-bit images depending on the maximum value of the image. -The filter is applied up to the image border, the neighboorhood used is +The filter is applied up to the image border, the neighborhood used is adjusted accordingly. The user may provide a mask image (same size as input image) where non zero values are the part of the image participating in the histogram computation. By default the entire image is filtered. diff --git a/skimage/filters/rank/tests/test_rank.py b/skimage/filters/rank/tests/test_rank.py index 0fda186c..78a97e56 100644 --- a/skimage/filters/rank/tests/test_rank.py +++ b/skimage/filters/rank/tests/test_rank.py @@ -87,7 +87,6 @@ def check_all(): def test_random_sizes(): # make sure the size is not a problem - niter = 10 elem = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]], dtype=np.uint8) for m, n in np.random.random_integers(1, 100, size=(10, 2)): mask = np.ones((m, n), dtype=np.uint8) From 1847b18fbba0998a585050d722d9a2167240c302 Mon Sep 17 00:00:00 2001 From: Juan Nunez-Iglesias Date: Thu, 22 Jan 2015 18:30:10 +1100 Subject: [PATCH 36/36] Test symmetric selems for rank/morpho comparison The handling of eccentric selems is different between these two modules, so the results are not expected to agree. (Why they agreed before this change is beyond me.) --- skimage/filters/rank/tests/test_rank.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/filters/rank/tests/test_rank.py b/skimage/filters/rank/tests/test_rank.py index 78a97e56..48628432 100644 --- a/skimage/filters/rank/tests/test_rank.py +++ b/skimage/filters/rank/tests/test_rank.py @@ -124,7 +124,7 @@ def test_compare_with_grey_dilation(): out = np.empty_like(image) mask = np.ones(image.shape, dtype=np.uint8) - for r in range(1, 20, 1): + for r in range(3, 20, 2): elem = np.ones((r, r), dtype=np.uint8) rank.maximum(image=image, selem=elem, out=out, mask=mask) cm = grey.dilation(image=image, selem=elem) @@ -138,7 +138,7 @@ def test_compare_with_grey_erosion(): out = np.empty_like(image) mask = np.ones(image.shape, dtype=np.uint8) - for r in range(1, 20, 1): + for r in range(3, 20, 2): elem = np.ones((r, r), dtype=np.uint8) rank.minimum(image=image, selem=elem, out=out, mask=mask) cm = grey.erosion(image=image, selem=elem)