wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-06-28 23:57:58 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1310 from jni/morpho-update

Browse Source 
Wrap ndimage morphology
This commit is contained in:
Stefan van der Walt
2015-01-22 00:06:23 -08:00
parent 9b16b35d96 1847b18fbb
commit d0366d2351
14 changed files with 601 additions and 244 deletions
Download Patch File Download Diff File Expand all files Collapse all files
bento.info
-3
View File
@@ -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
skimage/filters/rank/README.rst
+2 -2
View File
@@ -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.
@@ -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.
skimage/filters/rank/tests/test_rank.py
+7 -8
View File
@@ -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
@@ -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)
@@ -118,31 +117,31 @@ 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)
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 = 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)
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 = cmorph._erode(image=image, selem=elem)
cm = grey.erosion(image=image, selem=elem)
assert_equal(out, cm)
skimage/morphology/binary.py
+18 -41
View File
@@ -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 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.
@@ -35,27 +33,17 @@ 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]``.
"""
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.
@@ -81,26 +69,16 @@ 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]``.
"""
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 +112,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 +141,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
skimage/morphology/grey.py
+177 -61
View File
@@ -1,17 +1,134 @@
"""
Grayscale morphological operations
"""
from skimage import img_as_ubyte
from .misc import default_fallback
from . import cmorph
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']
@default_fallback
def _shift_selem(selem, shift_x, shift_y):
"""Shift the binary image `selem` in the left and/or up.
This only affects 2D 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
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, same shape and type as `selem`
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
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. Its parameters must
include at least `image`, `selem`, and `out`.
Returns
-------
func_out : callable
The same function, but correctly padding the input image before
applying the input 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.
@@ -24,7 +141,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,14 +152,14 @@ 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
-----
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
--------
@@ -62,16 +179,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.
@@ -96,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
@@ -123,17 +239,22 @@ 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)
# 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)
nd.grey_dilation(image, footprint=selem, output=out)
return out
@default_fallback
@default_selem
@pad_for_eccentric_selems
def opening(image, selem=None, out=None):
"""Return greyscale morphological opening of an image.
@@ -147,7 +268,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 +276,7 @@ def opening(image, selem=None, out=None):
Returns
-------
opening : uint8 array
opening : array, same shape and type as `image`
The result of the morphological opening.
Examples
@@ -176,17 +297,14 @@ 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
@pad_for_eccentric_selems
def closing(image, selem=None, out=None):
"""Return greyscale morphological closing of an image.
@@ -200,7 +318,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 +326,7 @@ def closing(image, selem=None, out=None):
Returns
-------
closing : uint8 array
closing : array, same shape and type as `image`
The result of the morphological closing.
Examples
@@ -229,17 +347,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 +366,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 +374,7 @@ def white_tophat(image, selem=None, out=None):
Returns
-------
opening : uint8 array
out : array, same shape and type as `image`
The result of the morphological white top hat.
Examples
@@ -281,16 +395,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.
@@ -312,7 +428,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
@@ -333,10 +449,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
skimage/morphology/misc.py
+6 -20
View File
@@ -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
skimage/morphology/setup.py
-3
View File
@@ -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'],
skimage/morphology/tests/test_binary.py
+8 -15
View File
@@ -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)
skimage/morphology/tests/test_grey.py
+39 -39
View File
@@ -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():
skimage/segmentation/boundaries.py
+198 -21
View File
@@ -1,41 +1,218 @@
from __future__ import division
import numpy as np
from ..morphology import dilation, square
from ..util import img_as_float
from scipy import ndimage as nd
from ..morphology import dilation, erosion, square
from ..util import img_as_float, view_as_windows, pad
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_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.
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.
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. 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, 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],
[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)
>>> 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)
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)
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:
boundaries = _find_boundaries_subpixel(label_img)
return boundaries
def mark_boundaries(image, label_img, color=(1, 1, 0),
outline_color=(0, 0, 0)):
outline_color=None, mode='outer', background_label=0):
"""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.
"""
if image.ndim == 2:
image = gray2rgb(image)
image = img_as_float(image, force_copy=True)
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``).
boundaries = find_boundaries(label_img)
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``.
"""
marked = img_as_float(image, force_copy=True)
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,
background=background_label)
if outline_color is not None:
outer_boundaries = dilation(boundaries.astype(np.uint8), square(2))
image[outer_boundaries != 0, :] = np.array(outline_color)
image[boundaries, :] = np.array(color)
return image
outlines = dilation(boundaries, square(3))
marked[outlines] = outline_color
marked[boundaries] = color
return marked
skimage/segmentation/tests/test_boundaries.py
+57 -27
View File
@@ -1,19 +1,22 @@
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))
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],
[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]])
@@ -24,36 +27,63 @@ 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],
[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)
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, 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)
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)
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 = 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()
skimage/util/__init__.py
+2 -1
View File
@@ -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']
skimage/util/arraypad.py
+41 -1
View File
@@ -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 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 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
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, copy=False)
crops = _validate_lengths(ar, crop_width)
slices = [slice(a, ar.shape[i] - b) for i, (a, b) in enumerate(crops)]
if copy:
cropped = np.array(ar[slices], order=order, copy=True)
else:
cropped = ar[slices]
return cropped
skimage/util/tests/test_arraypad.py
+46 -2
View File
@@ -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,49 @@ 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.flags.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)
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()