diff --git a/skimage/morphology/__init__.py b/skimage/morphology/__init__.py index efa2077d..d4c775eb 100644 --- a/skimage/morphology/__init__.py +++ b/skimage/morphology/__init__.py @@ -1,3 +1,5 @@ +from .binary import (binary_erosion, binary_dilation, binary_opening, + binary_closing) from .grey import * from .selem import * from .ccomp import label diff --git a/skimage/morphology/binary.py b/skimage/morphology/binary.py new file mode 100644 index 00000000..eef31a3d --- /dev/null +++ b/skimage/morphology/binary.py @@ -0,0 +1,133 @@ +import numpy as np +from scipy import ndimage + + +def binary_erosion(image, selem, out=None): + """Return fast binary morphological erosion of an image. + + This function returns the same result as greyscale erosion but performs + faster for binary images. + + Morphological erosion sets a pixel at (i,j) to the minimum over all pixels + in the neighborhood centered at (i,j). Erosion shrinks bright regions and + enlarges dark regions. + + Parameters + ---------- + image : ndarray + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + The array to store the result of the morphology. If None is + passed, a new array will be allocated. + + Returns + ------- + eroded : bool array + The result of the morphological erosion. + + """ + + out = ndimage.convolve(image > 0, selem, output=out, + mode='constant', cval=1) + return np.equal(out, np.sum(selem), out=out) + + +def binary_dilation(image, selem, out=None): + """Return fast binary morphological dilation of an image. + + This function returns the same result as greyscale dilation but performs + faster for binary images. + + Morphological dilation sets a pixel at (i,j) to the maximum over all pixels + in the neighborhood centered at (i,j). Dilation enlarges bright regions + and shrinks dark regions. + + Parameters + ---------- + + image : ndarray + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + The array to store the result of the morphology. If None, is + passed, a new array will be allocated. + + Returns + ------- + dilated : bool array + The result of the morphological dilation. + + """ + + out = ndimage.convolve(image > 0, selem, output=out, + mode='constant', cval=1) + return np.not_equal(out, 0, out=out) + + +def binary_opening(image, selem, out=None): + """Return fast binary morphological opening of an image. + + This function returns the same result as greyscale opening but performs + faster for binary images. + + The morphological opening on an image is defined as an erosion followed by + a dilation. Opening can remove small bright spots (i.e. "salt") and connect + small dark cracks. This tends to "open" up (dark) gaps between (bright) + features. + + Parameters + ---------- + image : ndarray + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + The array to store the result of the morphology. If None + is passed, a new array will be allocated. + + Returns + ------- + opening : bool array + The result of the morphological opening. + + """ + + eroded = binary_erosion(image, selem) + out = binary_dilation(eroded, selem, out=out) + return out + + +def binary_closing(image, selem, out=None): + """Return fast binary morphological closing of an image. + + This function returns the same result as greyscale closing but performs + faster for binary images. + + The morphological closing on an image is defined as a dilation followed by + an erosion. Closing can remove small dark spots (i.e. "pepper") and connect + small bright cracks. This tends to "close" up (dark) gaps between (bright) + features. + + Parameters + ---------- + image : ndarray + Image array. + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + The array to store the result of the morphology. If None, + is passed, a new array will be allocated. + + Returns + ------- + closing : bool array + The result of the morphological closing. + + """ + + dilated = binary_dilation(image, selem) + out = binary_erosion(dilated, selem, out=out) + return out diff --git a/skimage/morphology/cmorph.pyx b/skimage/morphology/cmorph.pyx index 6870b500..fb1e3af1 100644 --- a/skimage/morphology/cmorph.pyx +++ b/skimage/morphology/cmorph.pyx @@ -1,116 +1,118 @@ -""" -:author: Damian Eads, 2009 -:license: modified BSD -""" +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False -from __future__ import division import numpy as np - cimport numpy as np -cimport cython -from cpython cimport bool +from libc.stdlib cimport malloc, free -STREL_DTYPE = np.uint8 -ctypedef np.uint8_t STREL_DTYPE_t -IMAGE_DTYPE = np.uint8 -ctypedef np.uint8_t IMAGE_DTYPE_t +def dilate(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): -cdef inline int int_max(int a, int b): return a if a >= b else b -cdef inline int int_min(int a, int b): return a if a <= b else b + cdef int rows = image.shape[0] + cdef int cols = image.shape[1] + cdef int srows = selem.shape[0] + cdef int scols = selem.shape[1] -@cython.boundscheck(False) -def dilate(np.ndarray[IMAGE_DTYPE_t, ndim=2] image not None, - np.ndarray[IMAGE_DTYPE_t, ndim=2] selem not None, - np.ndarray[IMAGE_DTYPE_t, ndim=2] out, - bool shift_x, bool shift_y): - cdef int hw = selem.shape[0] // 2 - cdef int hh = selem.shape[1] // 2 - if shift_x: - hh -= 1 - if shift_y: - hw -= 1 + cdef int centre_r = int(selem.shape[0] / 2) - shift_y + cdef int centre_c = int(selem.shape[1] / 2) - shift_x - cdef int width = image.shape[0], height = image.shape[1] + image = np.ascontiguousarray(image) if out is None: - out = np.zeros([width, height], dtype=IMAGE_DTYPE) + out = np.zeros((rows, cols), dtype=np.uint8) + else: + out = np.ascontiguousarray(out) - assert out.shape[0] == image.shape[0] - assert out.shape[1] == image.shape[1] + cdef np.uint8_t* out_data = out.data + cdef np.uint8_t* image_data = image.data - cdef int x, y, ix, iy, cx, cy - cdef IMAGE_DTYPE_t max_so_far + cdef int r, c, rr, cc, s, value, local_max - cdef int sw = selem.shape[0], sh = selem.shape[1] + cdef int selem_num = np.sum(selem != 0) + cdef int* sr = malloc(selem_num * sizeof(int)) + cdef int* sc = malloc(selem_num * sizeof(int)) - cdef np.ndarray[np.int_t, ndim=2] xinc = np.zeros([sw, sh], dtype=np.int) - cdef np.ndarray[np.int_t, ndim=2] yinc = np.zeros([sw, sh], dtype=np.int) + s = 0 + for r in range(srows): + for c in range(scols): + if selem[r, c] != 0: + sr[s] = r - centre_r + sc[s] = c - centre_c + s += 1 - for x in range(sw): - for y in range(sh): - xinc[x, y] = (x - hw) - yinc[x, y] = (y - hh) + for r in range(rows): + for c in range(cols): + local_max = 0 + for s in range(selem_num): + rr = r + sr[s] + cc = c + sc[s] + if 0 <= rr < rows and 0 <= cc < cols: + value = image_data[rr * rows + cc] + if value > local_max: + local_max = value + out_data[r * cols + c] = local_max - for x in range(width): - for y in range(height): - max_so_far = 0 - for cx in range(0, sw): - for cy in range(0, sh): - ix = x + xinc[cx,cy] - iy = y + yinc[cx,cy] - if ix>=0 and iy>=0 and ix < width and iy < height \ - and selem[cx, cy] == 1 \ - and image[ix,iy] > max_so_far: - max_so_far = image[ix,iy] - out[x,y] = max_so_far + free(sr) + free(sc) return out -@cython.boundscheck(False) -def erode(np.ndarray[IMAGE_DTYPE_t, ndim=2] image not None, - np.ndarray[IMAGE_DTYPE_t, ndim=2] selem not None, - np.ndarray[IMAGE_DTYPE_t, ndim=2] out, - bool shift_x, bool shift_y): - cdef int hw = selem.shape[0] // 2 - cdef int hh = selem.shape[1] // 2 - if shift_x: - hh -= 1 - if shift_y: - hw -= 1 +def erode(np.ndarray[np.uint8_t, ndim=2] image, + np.ndarray[np.uint8_t, ndim=2] selem, + np.ndarray[np.uint8_t, ndim=2] out=None, + char shift_x=0, char shift_y=0): - cdef int width = image.shape[0], height = image.shape[1] + cdef int rows = image.shape[0] + cdef int cols = image.shape[1] + cdef int srows = selem.shape[0] + cdef int scols = selem.shape[1] + + cdef int centre_r = int(selem.shape[0] / 2) - shift_y + cdef int centre_c = int(selem.shape[1] / 2) - shift_x + + image = np.ascontiguousarray(image) if out is None: - out = np.zeros([width, height], dtype=IMAGE_DTYPE) + out = np.zeros((rows, cols), dtype=np.uint8) + else: + out = np.ascontiguousarray(out) - assert out.shape[0] == image.shape[0] - assert out.shape[1] == image.shape[1] + cdef np.uint8_t* out_data = out.data + cdef np.uint8_t* image_data = image.data - cdef int x, y, ix, iy, cx, cy - cdef IMAGE_DTYPE_t min_so_far + cdef int r, c, rr, cc, s, value, local_max - cdef int sw = selem.shape[0], sh = selem.shape[1] + cdef int selem_num = np.sum(selem != 0) + cdef int* sr = malloc(selem_num * sizeof(int)) + cdef int* sc = malloc(selem_num * sizeof(int)) - cdef np.ndarray[np.int_t, ndim=2] xinc = np.zeros([sw, sh], dtype=np.int) - cdef np.ndarray[np.int_t, ndim=2] yinc = np.zeros([sw, sh], dtype=np.int) + s = 0 + for r in range(srows): + for c in range(scols): + if selem[r, c] != 0: + sr[s] = r - centre_r + sc[s] = c - centre_c + s += 1 - for x in range(sw): - for y in range(sh): - xinc[x, y] = (x - hw) - yinc[x, y] = (y - hh) + for r in range(rows): + for c in range(cols): + local_min = 255 + for s in range(selem_num): + rr = r + sr[s] + cc = c + sc[s] + if 0 <= rr < rows and 0 <= cc < cols: + value = image_data[rr * rows + cc] + if value < local_min: + local_min = value - for x in range(width): - for y in range(height): - min_so_far = 255 - for cx in range(0, sw): - for cy in range(0, sh): - ix = x + xinc[cx,cy] - iy = y + yinc[cx,cy] - if ix>=0 and iy>=0 and ix < width \ - and iy < height and selem[cx, cy] == 1 \ - and image[ix,iy] < min_so_far: - min_so_far = image[ix,iy] - out[x,y] = min_so_far + out_data[r * cols + c] = local_min + + free(sr) + free(sc) return out diff --git a/skimage/morphology/grey.py b/skimage/morphology/grey.py index 309fb5be..5ef6b9af 100644 --- a/skimage/morphology/grey.py +++ b/skimage/morphology/grey.py @@ -6,11 +6,11 @@ __docformat__ = 'restructuredtext en' import warnings - import numpy as np - import skimage +from . import cmorph + __all__ = ['erosion', 'dilation', 'opening', 'closing', 'white_tophat', 'black_tophat', 'greyscale_erode', 'greyscale_dilate', @@ -28,15 +28,12 @@ def erosion(image, selem, out=None, shift_x=False, shift_y=False): Parameters ---------- image : ndarray - Image array. - + Image array. selem : ndarray - The neighborhood expressed as a 2-D array of 1's and 0's. - + The neighborhood expressed as a 2-D array of 1's and 0's. out : ndarray - The array to store the result of the morphology. If None is - passed, a new array will be allocated. - + The array to store the result of the morphology. If None is + passed, a new array will be allocated. shift_x, shift_y : bool shift structuring element about center point. This only affects eccentric structuring elements (i.e. selem with even numbered sides). @@ -44,7 +41,7 @@ def erosion(image, selem, out=None, shift_x=False, shift_y=False): Returns ------- eroded : uint8 array - The result of the morphological erosion. + The result of the morphological erosion. Examples -------- @@ -63,17 +60,13 @@ def erosion(image, selem, 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 = skimage.img_as_ubyte(image) - - try: - import skimage.morphology.cmorph as cmorph - out = cmorph.erode(image, selem, out=out, - shift_x=shift_x, shift_y=shift_y) - return out - except ImportError: - raise ImportError("cmorph extension not available.") + selem = skimage.img_as_ubyte(selem) + return cmorph.erode(image, selem, out=out, + shift_x=shift_x, shift_y=shift_y) def dilation(image, selem, out=None, shift_x=False, shift_y=False): @@ -87,15 +80,12 @@ def dilation(image, selem, out=None, shift_x=False, shift_y=False): ---------- image : ndarray - Image array. - + Image array. selem : ndarray - The neighborhood expressed as a 2-D array of 1's and 0's. - + The neighborhood expressed as a 2-D array of 1's and 0's. out : ndarray - The array to store the result of the morphology. If None, is - passed, a new array will be allocated. - + The array to store the result of the morphology. If None, is + passed, a new array will be allocated. shift_x, shift_y : bool shift structuring element about center point. This only affects eccentric structuring elements (i.e. selem with even numbered sides). @@ -103,7 +93,7 @@ def dilation(image, selem, out=None, shift_x=False, shift_y=False): Returns ------- dilated : uint8 array - The result of the morphological dilation. + The result of the morphological dilation. Examples -------- @@ -122,17 +112,13 @@ def dilation(image, selem, 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 = skimage.img_as_ubyte(image) - - try: - from . import cmorph - out = cmorph.dilate(image, selem, out=out, - shift_x=shift_x, shift_y=shift_y) - return out - except ImportError: - raise ImportError("cmorph extension not available.") + selem = skimage.img_as_ubyte(selem) + return cmorph.dilate(image, selem, out=out, + shift_x=shift_x, shift_y=shift_y) def opening(image, selem, out=None): @@ -146,19 +132,17 @@ def opening(image, selem, out=None): Parameters ---------- image : ndarray - Image array. - + Image array. selem : ndarray - The neighborhood expressed as a 2-D array of 1's and 0's. - + The neighborhood expressed as a 2-D array of 1's and 0's. out : ndarray - The array to store the result of the morphology. If None - is passed, a new array will be allocated. + The array to store the result of the morphology. If None + is passed, a new array will be allocated. Returns ------- opening : uint8 array - The result of the morphological opening. + The result of the morphological opening. Examples -------- @@ -177,6 +161,7 @@ def opening(image, selem, 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 @@ -197,19 +182,17 @@ def closing(image, selem, out=None): Parameters ---------- image : ndarray - Image array. - + Image array. selem : ndarray - The neighborhood expressed as a 2-D array of 1's and 0's. - + The neighborhood expressed as a 2-D array of 1's and 0's. out : ndarray - The array to store the result of the morphology. If None, - is passed, a new array will be allocated. + The array to store the result of the morphology. If None, + is passed, a new array will be allocated. Returns ------- closing : uint8 array - The result of the morphological closing. + The result of the morphological closing. Examples -------- @@ -228,6 +211,7 @@ def closing(image, selem, 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 @@ -247,19 +231,17 @@ def white_tophat(image, selem, out=None): Parameters ---------- image : ndarray - Image array. - + Image array. selem : ndarray - The neighborhood expressed as a 2-D array of 1's and 0's. - + The neighborhood expressed as a 2-D array of 1's and 0's. out : ndarray - The array to store the result of the morphology. If None - is passed, a new array will be allocated. + The array to store the result of the morphology. If None + is passed, a new array will be allocated. Returns ------- opening : uint8 array - The result of the morphological white top hat. + The result of the morphological white top hat. Examples -------- @@ -298,14 +280,12 @@ def black_tophat(image, selem, out=None): Parameters ---------- image : ndarray - Image array. - + Image array. selem : ndarray - The neighborhood expressed as a 2-D array of 1's and 0's. - + The neighborhood expressed as a 2-D array of 1's and 0's. out : ndarray - The array to store the result of the morphology. If None - is passed, a new array will be allocated. + The array to store the result of the morphology. If None + is passed, a new array will be allocated. Returns ------- @@ -329,6 +309,7 @@ def black_tophat(image, selem, out=None): [0, 0, 0, 0, 0]], dtype='uint8') """ + if image is out: raise NotImplementedError("Cannot perform white top hat in place.") image = skimage.img_as_ubyte(image) diff --git a/skimage/morphology/tests/test_grey.py b/skimage/morphology/tests/test_grey.py index d7242e8f..6003425f 100644 --- a/skimage/morphology/tests/test_grey.py +++ b/skimage/morphology/tests/test_grey.py @@ -5,11 +5,11 @@ from numpy import testing import skimage from skimage import data_dir -from skimage.morphology import grey -from skimage.morphology import selem +from skimage.morphology import binary, grey, selem lena = np.load(os.path.join(data_dir, 'lena_GRAY_U8.npy')) +bw_lena = lena > 0.4 class TestMorphology(): @@ -154,5 +154,33 @@ class TestDTypes(): self._test_image(image) +def test_binary_erosion(): + strel = selem.square(3) + binary_res = binary.binary_erosion(bw_lena, strel) + grey_res = grey.erosion(bw_lena, strel) + testing.assert_array_equal(binary_res, grey_res) + + +def test_binary_dilation(): + strel = selem.square(3) + binary_res = binary.binary_dilation(bw_lena, strel) + grey_res = grey.dilation(bw_lena, strel) + testing.assert_array_equal(binary_res, grey_res) + + +def test_binary_closing(): + strel = selem.square(3) + binary_res = binary.binary_closing(bw_lena, strel) + grey_res = grey.closing(bw_lena, strel) + testing.assert_array_equal(binary_res, grey_res) + + +def test_binary_opening(): + strel = selem.square(3) + binary_res = binary.binary_opening(bw_lena, strel) + grey_res = grey.opening(bw_lena, strel) + testing.assert_array_equal(binary_res, grey_res) + + if __name__ == '__main__': testing.run_module_suite()