Finally, I've checked in my additions (hopefully) the right way.

This commit is contained in:
Damian Eads
2009-10-27 01:13:54 -07:00
parent 3bf33a0a46
commit 382fa9c788
5 changed files with 4496 additions and 0 deletions
+2
View File
@@ -0,0 +1,2 @@
from grey import *
from selem import *
File diff suppressed because it is too large Load Diff
+99
View File
@@ -0,0 +1,99 @@
"""
:author: Damian Eads, 2009
:license: modified BSD
"""
from __future__ import division
import numpy as np
cimport numpy as np
cimport cython
STREL_DTYPE = np.uint8
ctypedef np.uint8_t STREL_DTYPE_t
IMAGE_DTYPE = np.uint8
ctypedef np.uint8_t IMAGE_DTYPE_t
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
@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):
cdef int hw = selem.shape[0] / 2
cdef int hh = selem.shape[1] / 2
cdef int width = image.shape[0], height = image.shape[1]
if out is None:
out = np.zeros([width, height], dtype=IMAGE_DTYPE)
assert out.shape[0] == image.shape[0]
assert out.shape[1] == image.shape[1]
cdef int x, y, ix, iy, cx, cy
cdef IMAGE_DTYPE_t max_so_far
cdef int sw = selem.shape[0], sh = selem.shape[1]
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)
for x in range(sw):
for y in range(sh):
xinc[x, y] = (x - hw)
yinc[x, y] = (y - hh)
for x in range(width):
for y in range(height):
max_so_far = 0
#for cx in range(int_max(0, x - hw), int_min(x + hw, out.shape[0])):
#for cy in range(int_max(0, y - hh), int_min(y + hh, out.shape[1])):
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
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):
cdef int hw = selem.shape[0] / 2
cdef int hh = selem.shape[1] / 2
cdef int width = image.shape[0], height = image.shape[1]
if out is None:
out = np.zeros([width, height], dtype=IMAGE_DTYPE)
assert out.shape[0] == image.shape[0]
assert out.shape[1] == image.shape[1]
cdef int x, y, ix, iy, cx, cy
cdef IMAGE_DTYPE_t min_so_far
cdef int sw = selem.shape[0], sh = selem.shape[1]
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)
for x in range(sw):
for y in range(sh):
xinc[x, y] = (x - hw)
yinc[x, y] = (y - hh)
for x in range(width):
for y in range(height):
min_so_far = 255
#for cx in range(int_max(0, x - hw), int_min(x + hw, out.shape[0])):
#for cy in range(int_max(0, y - hh), int_min(y + hh, out.shape[1])):
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
return out
+189
View File
@@ -0,0 +1,189 @@
"""
:author: Damian Eads, 2009
:license: modified BSD
"""
#__all__ = ['square', 'disk', 'diamond', 'line', 'ball', ]
__docformat__ = 'restructuredtext en'
import numpy as np
#from scipy.fftpack import fftshift, ifftshift
eps = np.finfo(float).eps
def greyscale_erode(image, selem, out=None):
"""
Performs a greyscale morphological erosion on an image given a
structuring element. The eroded pixel at (i,j) is the minimum
over all pixels in the neighborhood centered at (i,j).
Parameters
----------
image : ndarray
The image as an ndarray.
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 : ndarray
The result of the morphological erosion.
"""
if image is out:
raise NotImplementedError("In-place morphological erosion not supported!")
try:
import cmorph
out = cmorph.erode(image, selem, out=out)
return out;
except ImportError:
raise ImportError("cmorph extension not available.")
def greyscale_dilate(image, selem, out=None):
"""
Performs a greyscale morphological dilation on an image given a
structuring element. The dilated pixel at (i,j) is the maximum
over all pixels in the neighborhood centered at (i,j).
Parameters
----------
image : ndarray
The image as an ndarray.
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 : ndarray
The result of the morphological dilation.
"""
if image is out:
raise NotImplementedError("In-place morphological dilation not supported!")
try:
import cmorph
out = cmorph.dilate(image, selem, out=out)
return out;
except ImportError:
raise ImportError("cmorph extension not available.")
def greyscale_open(image, selem, out=None):
"""
Performs a greyscale morphological opening on an image given a
structuring element defined as a erosion followed by a dilation.
Parameters
----------
image : ndarray
The image as an ndarray.
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 : ndarray
The result of the morphological opening.
"""
eroded = greyscale_erode(image, selem)
out = greyscale_dilate(eroded, selem, out=out)
return out
def greyscale_close(image, selem, out=None):
"""
Performs a greyscale morphological closing on an image given a
structuring element defined as a dilation followed by an erosion.
Parameters
----------
image : ndarray
The image as an ndarray.
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 : ndarray
The result of the morphological opening.
"""
dilated = greyscale_dilate(image, selem)
out = greyscale_erode(dilated, selem, out=out)
return out
def greyscale_white_top_hat(image, selem, out=None):
"""
Applies a white top hat on an image given a structuring element.
Parameters
----------
image : ndarray
The image as an ndarray.
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 : ndarray
The result of the morphological white top hat.
"""
if image is out:
raise NotImplementedError("Cannot perform white top hat in place.")
eroded = greyscale_erode(image, selem)
out = greyscale_dilate(eroded, selem, out=out)
out = image - out
return out
def greyscale_black_top_hat(image, selem, out=None):
"""
Applies a black top hat on an image given a structuring element.
Parameters
----------
image : ndarray
The image as an ndarray.
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 : ndarray
The result of the black top filter.
"""
if image is out:
raise NotImplementedError("Cannot perform white top hat in place.")
dilated = greyscale_dilate(image, selem)
out = greyscale_erode(dilated, selem, out=out)
out = out - image
if image is out:
raise NotImplementedError("Cannot perform black top hat in place.")
return out
+209
View File
@@ -0,0 +1,209 @@
"""
:author: Damian Eads, 2009
:license: modified BSD
"""
import numpy as np
def square(width, dtype='uint8'):
"""
Generates a flat, square-shaped structuring element. Every pixel
along the perimeter has a chessboard distance no greater than radius
(radius=floor(width/2)) pixels.
Parameters
----------
width : int
The width and height of the square
Additional Parameters
---------------------
dtype : string
The data type of the structuring element.
Returns
-------
selem : ndarray
A structuring element consisting only of ones, i.e. every
pixel belongs to the neighborhood.
"""
return np.ones((width, width), dtype=dtype)
def rectangle(width, height, dtype='uint8'):
"""
Generates a flat, rectangular-shaped structuring element of a
given width and height. Every pixel in the rectangle belongs
to the neighboorhood.
Parameters
----------
width : int
The width of the rectangle
height : int
The height of the rectangle
Additional Parameters
---------------------
dtype : string
The data type of the structuring element.
Returns
-------
selem : ndarray
A structuring element consisting only of ones, i.e. every
pixel belongs to the neighborhood.
"""
return np.ones((width, height), dtype=dtype)
def diamond(radius, dtype='uint8'):
"""
Generates a flat, diamond-shaped structuring element of a given
radius. A pixel is part of the neighborhood (i.e. labeled 1) iff
the city block/manhattan distance between it and the center of the
neighborhood is no greater than radius.
*Parameters*:
radius : string
The radius of the disk-shaped structuring element.
dtype : string
The data type of the structuring element.
*Returns*:
selem : ndarray
The structuring element where elements of the neighborhood
are 1 and 0 otherwise.
"""
half = radius
(I, J) = np.meshgrid(xrange(0, radius*2+1), xrange(0, radius*2+1))
s = np.abs(I-half)+np.abs(J-half)
return np.array(s <= radius, dtype=dtype)
def disk(radius, N=0, dtype='uint8'):
"""
Generates a flat, disk-shaped structuring element of a given radius.
A pixel is within the neighborhood iff the euclidean distance between
it and the origin is no greater than a radius.
Parameters
----------
radius : string
The radius of the disk-shaped structuring element.
dtype : string
The data type of the structuring element.
Returns
-------
selem : ndarray
The structuring element where elements of the neighborhood
are 1 and 0 otherwise.
"""
half = radius
(I, J) = np.meshgrid(xrange(0, radius*2+1), xrange(0, radius*2+1))
if N == 0:
s = (I-half)**2.+(J-half)**2.
#print s
else:
raise NotImplementedError("scikits.image.morphology.disk: approximations not implemented. Try N=0 for now.")
return np.array(s <= radius * radius, dtype=dtype)
def ellipse(size, angle, ratio=0.5, dtype='uint8'):
"""
Generates an elliptically-shaped structuring element of a given
angle, ratio, and kernel size.
Parameters
----------
size : int
The half-width of the kernel. The kernel size is 2*size+1.
angle : float
The angle of rotation of the ellipse in radians.
ratio : float
The aspect ratio of the ellipse.
dtype : string
The data type of the structuring element.
Returns
-------
selem : ndarray
The structuring element where elements of the neighborhood
are 1 and 0 otherwise.
"""
structure = np.zeros((2*size+1, 2*size+1), dtype=dtype)
a = np.matrix([np.cos(angle), np.sin(angle)])
b = np.matrix([-np.sin(angle), np.cos(angle)])
aspect = a.T * a + b.T * b / ratio
for y in xrange(-size, size+ 1):
for x in xrange(-size, size + 1):
i = x+size
j = y+size
v = np.matrix([x,y], dtype='f') / float(size)
n = v * aspect * v.T
if n < 1:
structure[i, j] = 1
return structure
def strel(shape='disk', N=0, radius=3, width=3, height=3, angle=0., length=3, dtype='uint8', out=None):
"""
Generates a structuring element for greyscale or binary morphology.
The interface of this function is similar to MATLAB(TM)'s strel function.
Parameters
----------
shape : string
A string identifier for the shape of the structuring element,
which can be any of the following: 'arbitrary', 'ball',
'diamond', 'disk', 'pair', 'rectangle', 'square'.
N : int
When non-zero, the number of lines to approximate the
structuring element. (not implemented)
radius : int
The radius for disk or diamond structuring elements.
width : int
The height for square, ball or rectangle-shaped structuring elements.
height : int
The height for ball or rectangle-shaped structuring elements.
size : int
The half-width of an elliptical structuring element.
aspect : float
The aspect ratio of an ellipse.
Returns
-------
neighborhood : ndarray
The structuring element.
"""
shape = shape.lower().strip()
if shape == 'disk':
return disk(radius, dtype=dtype)
elif shape == 'diamond':
return diamond(radius, dtype=dtype)
elif shape == 'square':
return square(width, dtype=dtype)
elif shape == 'rectangle':
return rectangle(width, height, dtype=dtype)
elif shape == 'ellipse':
return ellipse(size, angle, ratio, dtype=dtype)
else:
raise ValueError("Unknown structuring element type '%s'" % shape)