Merge pull request #287 from ahojnnes/morph

ENH: Improve performance of erosion and dilation.
This commit is contained in:
Stefan van der Walt
2012-09-02 04:36:19 -07:00
5 changed files with 295 additions and 149 deletions
+2
View File
@@ -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
+133
View File
@@ -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
+88 -86
View File
@@ -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 = <np.uint8_t*>out.data
cdef np.uint8_t* image_data = <np.uint8_t*>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 = <int*>malloc(selem_num * sizeof(int))
cdef int* sc = <int*>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 = <np.uint8_t*>out.data
cdef np.uint8_t* image_data = <np.uint8_t*>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 = <int*>malloc(selem_num * sizeof(int))
cdef int* sc = <int*>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
+42 -61
View File
@@ -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)
+30 -2
View File
@@ -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()