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

Merge pull request #215 from tonysyu/morph-reconstruction

Browse Source 
ENH: Add morphological reconstruction.
This commit is contained in:
Stefan van der Walt
2012-08-25 07:42:36 -07:00
parent 3c14898fda 666d6f4f09
commit 544fec7560
8 changed files with 560 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files
DEPENDS.txt
+1 -1
View File
@@ -2,7 +2,7 @@ Build Requirements
------------------
* `Python >= 2.5 <http://python.org>`__
* `Numpy >= 1.6 <http://numpy.scipy.org/>`__
* `Cython >= 0.15 <http://www.cython.org/>`__
* `Cython >= 0.16 <http://www.cython.org/>`__
`Matplotlib >= 1.0 <http://matplotlib.sf.net>`__ is needed to generate the
examples in the documentation.
doc/examples/plot_holes_and_peaks.py
+86
View File
@@ -0,0 +1,86 @@
"""
===============================
Filling holes and finding peaks
===============================
In this example, we fill holes (i.e. isolated, dark spots) in an image using
morphological reconstruction by erosion. Erosion expands the minimal values of
the seed image until it encounters a mask image. Thus, the seed image and mask
image represent the maximum and minimum possible values of the reconstructed
image.
We start with an image containing both peaks and holes:
"""
import matplotlib.pyplot as plt
from skimage import data
from skimage.exposure import rescale_intensity
image = data.moon()
# Rescale image intensity so that we can see dim features.
image = rescale_intensity(image, in_range=(50, 200))
# convenience function for plotting images
def imshow(image, **kwargs):
plt.figure(figsize=(5, 4))
plt.imshow(image, **kwargs)
plt.axis('off')
imshow(image)
plt.title('original image')
"""
.. image:: PLOT2RST.current_figure
Now we need to create the seed image, where the minima represent the starting
points for erosion. To fill holes, we initialize the seed image to the maximum
value of the original image. Along the borders, however, we use the original
values of the image. These border pixels will be the starting points for the
erosion process. We then limit the erosion by setting the mask to the values
of the original image.
"""
import numpy as np
from skimage.morphology import reconstruction
seed = np.copy(image)
seed[1:-1, 1:-1] = image.max()
mask = image
filled = reconstruction(seed, mask, method='erosion')
imshow(filled, vmin=image.min(), vmax=image.max())
plt.title('after filling holes')
"""
.. image:: PLOT2RST.current_figure
As shown above, eroding inward from the edges removes holes, since (by
definition) holes are surrounded by pixels of brighter value. Finally, we can
isolate the dark regions by subtracting the reconstructed image from the
original image.
"""
imshow(image - filled)
plt.title('holes')
"""
.. image:: PLOT2RST.current_figure
Alternatively, we can find bright spots in an image using morphological
reconstruction by dilation. Dilation is the inverse of erosion and expands the
*maximal* values of the seed image until it encounters a mask image. Since this
is an inverse operation, we initialize the seed image to the minimum image
intensity instead of the maximum. The remainder of the process is the same.
"""
seed = np.copy(image)
seed[1:-1, 1:-1] = image.min()
rec = reconstruction(seed, mask, method='dilation')
imshow(image - rec)
plt.title('peaks')
plt.show()
"""
.. image:: PLOT2RST.current_figure
"""
doc/examples/plot_regional_maxima.py
+113
View File
@@ -0,0 +1,113 @@
"""
=========================
Filtering regional maxima
=========================
Here, we use morphological reconstruction to create a background image, which
we can subtract from the original image to isolate bright features (regional
maxima).
First we try reconstruction by dilation starting at the edges of the image. We
initialize a seed image to the minimum intensity of the image, and set its
border to be the pixel values in the original image. These maximal pixels will
get dilated in order to reconstruct the background image.
"""
import numpy as np
from skimage import data
from skimage import img_as_float
from skimage.morphology import reconstruction
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt
# Convert to float: Important for subtraction later which won't work with uint8
image = img_as_float(data.coins())
image = gaussian_filter(image, 1)
seed = np.copy(image)
seed[1:-1, 1:-1] = image.min()
mask = image
dilated = reconstruction(seed, mask, method='dilation')
"""
Subtracting the dilated image leaves an image with just the coins and a flat,
black background, as shown below.
"""
fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 2.5))
ax1.imshow(image)
ax1.set_title('original image')
ax1.axis('off')
ax2.imshow(dilated, vmin=image.min(), vmax=image.max())
ax2.set_title('dilated')
ax2.axis('off')
ax3.imshow(image - dilated)
ax3.set_title('image - dilated')
ax3.axis('off')
plt.tight_layout()
"""
.. image:: PLOT2RST.current_figure
Although the features (i.e. the coins) are clearly isolated, the coins
surrounded by a bright background in the original image are dimmer in the
subtracted image. We can attempt to correct this using a different seed image.
Instead of creating a seed image with maxima along the image border, we can use
the features of the image itself to seed the reconstruction process. Here, the
seed image is the original image minus a fixed value, ``h``.
"""
h = 0.4
seed = image - h
dilated = reconstruction(seed, mask, method='dilation')
hdome = image - dilated
"""
To get a feel for the reconstruction process, we plot the intensity of the
mask, seed, and dilated images along a slice of the image (indicated by red
line).
"""
fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 2.5))
yslice = 197
ax1.plot(mask[yslice], '0.5', label='mask')
ax1.plot(seed[yslice], 'k', label='seed')
ax1.plot(dilated[yslice], 'r', label='dilated')
ax1.set_ylim(-0.2, 2)
ax1.set_title('image slice')
ax1.set_xticks([])
ax1.legend()
ax2.imshow(dilated, vmin=image.min(), vmax=image.max())
ax2.axhline(yslice, color='r', alpha=0.4)
ax2.set_title('dilated')
ax2.axis('off')
ax3.imshow(hdome)
ax3.axhline(yslice, color='r', alpha=0.4)
ax3.set_title('image - dilated')
ax3.axis('off')
plt.tight_layout()
plt.show()
"""
.. image:: PLOT2RST.current_figure
As you can see in the image slice, each coin is given a different baseline
intensity in the reconstructed image; this is because we used the local
intensity (shifted by ``h``) as a seed value. As a result, the coins in the
subtracted image have similar pixel intensities. The final result is known as
the h-dome of an image since this tends to isolate regional maxima of height
``h``. This operation is particularly useful when your images are unevenly
illuminated.
"""
skimage/morphology/__init__.py
+1
View File
@@ -4,3 +4,4 @@ from .ccomp import label
from .watershed import watershed, is_local_maximum
from ._skeletonize import skeletonize, medial_axis
from .convex_hull import convex_hull_image
from .greyreconstruct import reconstruction
skimage/morphology/_greyreconstruct.pyx
+81
View File
@@ -0,0 +1,81 @@
"""
`reconstruction_loop` originally part of CellProfiler, code licensed under both GPL and BSD licenses.
Website: http://www.cellprofiler.org
Copyright (c) 2003-2009 Massachusetts Institute of Technology
Copyright (c) 2009-2011 Broad Institute
All rights reserved.
Original author: Lee Kamentsky
"""
cimport cython
@cython.boundscheck(False)
def reconstruction_loop(unsigned int[:] ranks, int[:] prev, int[:] next,
int[:] strides, int current_idx, int image_stride):
"""The inner loop for reconstruction.
This algorithm uses the rank-order of pixels. If low intensity pixels have
a low rank and high intensity pixels have a high rank, then this loop
performs reconstruction by dilation. If this ranking is reversed, the
result is reconstruction by erosion.
For each pixel in the seed image, check its neighbors. If its neighbor's
rank is below that of the current pixel, replace the neighbor's rank with
the rank of the current pixel. This dilation is limited by the mask, i.e.
the rank at each pixel cannot exceed the mask as that pixel.
Parameters
----------
ranks : array
The rank order of the flattened seed and mask images.
prev, next: arrays
Indices of previous and next pixels in rank sorted order.
strides : array
Strides to neighbors of the current pixel.
current_idx : int
Index of highest-ranked pixel used as starting point in loop.
image_stride : int
Stride between seed image and mask image in `ranks`.
"""
cdef unsigned int neighbor_rank, current_rank, mask_rank
cdef int i, current_link, neighbor_idx, nprev, nnext
cdef int nstrides = strides.shape[0]
while current_idx != -1:
if current_idx < image_stride:
current_rank = ranks[current_idx]
if current_rank == 0:
break
for i in range(nstrides):
neighbor_idx = current_idx + strides[i]
neighbor_rank = ranks[neighbor_idx]
# Only propagate neighbors ranked below the current rank
if neighbor_rank < current_rank:
mask_rank = ranks[neighbor_idx + image_stride]
# Only propagate neighbors ranked below the mask rank
if neighbor_rank < mask_rank:
# Raise the neighbor to the mask rank if
# the mask ranked below the current rank
if mask_rank < current_rank:
current_link = neighbor_idx + image_stride
ranks[neighbor_idx] = mask_rank
else:
current_link = current_idx
ranks[neighbor_idx] = current_rank
# unlink the neighbor
nprev = prev[neighbor_idx]
nnext = next[neighbor_idx]
next[nprev] = nnext
if nnext != -1:
prev[nnext] = nprev
# link to the neighbor after the current link
nnext = next[current_link]
next[neighbor_idx] = nnext
prev[neighbor_idx] = current_link
if nnext >= 0:
prev[nnext] = neighbor_idx
next[current_link] = neighbor_idx
current_idx = next[current_idx]
skimage/morphology/greyreconstruct.py
+193
View File
@@ -0,0 +1,193 @@
"""
This morphological reconstruction routine was adapted from CellProfiler, code
licensed under both GPL and BSD licenses.
Website: http://www.cellprofiler.org
Copyright (c) 2003-2009 Massachusetts Institute of Technology
Copyright (c) 2009-2011 Broad Institute
All rights reserved.
Original author: Lee Kamentsky
"""
import numpy as np
from skimage.filter.rank_order import rank_order
def reconstruction(seed, mask, method='dilation', selem=None, offset=None):
"""Perform a morphological reconstruction of an image.
Morphological reconstruction by dilation is similar to basic morphological
dilation: high-intensity values will replace nearby low-intensity values.
The basic dilation operator, however, uses a structuring element to
determine how far a value in the input image can spread. In contrast,
reconstruction uses two images: a "seed" image, which specifies the values
that spread, and a "mask" image, which gives the maximum allowed value at
each pixel. The mask image, like the structuring element, limits the spread
of high-intensity values. Reconstruction by erosion is simply the inverse:
low-intensity values spread from the seed image and are limited by the mask
image, which represents the minimum allowed value.
Alternatively, you can think of reconstruction as a way to isolate the
connected regions of an image. For dilation, reconstruction connects
regions marked by local maxima in the seed image: neighboring pixels
less-than-or-equal-to those seeds are connected to the seeded region.
Local maxima with values larger than the seed image will get truncated to
the seed value.
Parameters
----------
seed : ndarray
The seed image (a.k.a. marker image), which specifies the values that
are dilated or eroded.
mask : ndarray
The maximum (dilation) / minimum (erosion) allowed value at each pixel.
method : {'dilation'|'erosion'}
Perform reconstruction by dilation or erosion. In dilation (or
erosion), the seed image is dilated (or eroded) until limited by the
mask image. For dilation, each seed value must be less than or equal
to the corresponding mask value; for erosion, the reverse is true.
selem : ndarray
The neighborhood expressed as a 2-D array of 1's and 0's.
Returns
-------
reconstructed : ndarray
The result of morphological reconstruction.
Examples
--------
>>> import numpy as np
>>> from skimage.morphology import reconstruction
First, we create a sinusoidal mask image w/ peaks at middle and ends.
>>> x = np.linspace(0, 4 * np.pi)
>>> y_mask = np.cos(x)
Then, we create a seed image initialized to the minimum mask value (for
reconstruction by dilation, min-intensity values don't spread) and add
"seeds" to the left and right peak, but at a fraction of peak value (1).
>>> y_seed = y_mask.min() * np.ones_like(x)
>>> y_seed[0] = 0.5
>>> y_seed[-1] = 0
>>> y_rec = reconstruction(y_seed, y_mask)
The reconstructed image (or curve, in this case) is exactly the same as the
mask image, except that the peaks are truncated to 0.5 and 0. The middle
peak disappears completely: Since there were no seed values in this peak
region, its reconstructed value is truncated to the surrounding value (-1).
As a more practical example, we try to extract the bright features of an
image by subtracting a background image created by reconstruction.
>>> y, x = np.mgrid[:20:0.5, :20:0.5]
>>> bumps = np.sin(x) + np.sin(y)
To create the background image, set the mask image to the original image,
and the seed image to the original image with an intensity offset, `h`.
>>> h = 0.3
>>> seed = bumps - h
>>> background = reconstruction(seed, bumps)
The resulting reconstructed image looks exactly like the original image,
but with the peaks of the bumps cut off. Subtracting this reconstructed
image from the original image leaves just the peaks of the bumps
>>> hdome = bumps - background
This operation is known as the h-dome of the image and leaves features
of height `h` in the subtracted image.
Notes
-----
The algorithm is taken from:
[1] Robinson, "Efficient morphological reconstruction: a downhill filter",
Pattern Recognition Letters 25 (2004) 1759-1767.
Applications for greyscale reconstruction are discussed in:
[2] Vincent, L., "Morphological Grayscale Reconstruction in Image Analysis:
Applications and Efficient Algorithms", IEEE Transactions on Image
Processing (1993)
[3] Soille, P., "Morphological Image Analysis: Principles and Applications",
Chapter 6, 2nd edition (2003), ISBN 3540429883.
"""
assert tuple(seed.shape) == tuple(mask.shape)
if method == 'dilation' and np.any(seed > mask):
raise ValueError("Intensity of seed image must be less than that "
"of the mask image for reconstruction by dilation.")
elif method == 'erosion' and np.any(seed < mask):
raise ValueError("Intensity of seed image must be greater than that "
"of the mask image for reconstruction by erosion.")
try:
from ._greyreconstruct import reconstruction_loop
except ImportError:
raise ImportError("_greyreconstruct extension not available.")
if selem is None:
selem = np.ones([3] * seed.ndim, dtype=bool)
else:
selem = selem.copy()
if offset == None:
if not all([d % 2 == 1 for d in selem.shape]):
ValueError("Footprint dimensions must all be odd")
offset = np.array([d / 2 for d in selem.shape])
# Cross out the center of the selem
selem[[slice(d, d + 1) for d in offset]] = False
# Make padding for edges of reconstructed image so we can ignore boundaries
padding = (np.array(selem.shape) / 2).astype(int)
dims = np.zeros(seed.ndim + 1, dtype=int)
dims[1:] = np.array(seed.shape) + 2 * padding
dims[0] = 2
inside_slices = [slice(p, -p) for p in padding]
# Set padded region to minimum image intensity and mask along first axis so
# we can interleave image and mask pixels when sorting.
if method == 'dilation':
pad_value = np.min(seed)
elif method == 'erosion':
pad_value = np.max(seed)
images = np.ones(dims) * pad_value
images[[0] + inside_slices] = seed
images[[1] + inside_slices] = mask
# Create a list of strides across the array to get the neighbors within
# a flattened array
value_stride = np.array(images.strides[1:]) / images.dtype.itemsize
image_stride = images.strides[0] / images.dtype.itemsize
selem_mgrid = np.mgrid[[slice(-o, d - o)
for d, o in zip(selem.shape, offset)]]
selem_offsets = selem_mgrid[:, selem].transpose()
nb_strides = np.array([np.sum(value_stride * selem_offset)
for selem_offset in selem_offsets], np.int32)
images = images.flatten()
# Erosion goes smallest to largest; dilation goes largest to smallest.
index_sorted = np.argsort(images).astype(np.int32)
if method == 'dilation':
index_sorted = index_sorted[::-1]
# Make a linked list of pixels sorted by value. -1 is the list terminator.
prev = -np.ones(len(images), np.int32)
next = -np.ones(len(images), np.int32)
prev[index_sorted[1:]] = index_sorted[:-1]
next[index_sorted[:-1]] = index_sorted[1:]
# Cython inner-loop compares the rank of pixel values.
if method == 'dilation':
value_rank, value_map = rank_order(images)
elif method == 'erosion':
value_rank, value_map = rank_order(-images)
value_map = -value_map
start = index_sorted[0]
reconstruction_loop(value_rank, prev, next, nb_strides, start, image_stride)
# Reshape reconstructed image to original image shape and remove padding.
rec_img = value_map[value_rank[:image_stride]]
rec_img.shape = np.array(seed.shape) + 2 * padding
return rec_img[inside_slices]
skimage/morphology/setup.py
+3
View File
@@ -18,6 +18,7 @@ def configuration(parent_package='', top_path=None):
cython(['_skeletonize_cy.pyx'], working_path=base_path)
cython(['_pnpoly.pyx'], working_path=base_path)
cython(['_convex_hull.pyx'], working_path=base_path)
cython(['_greyreconstruct.pyx'], working_path=base_path)
config.add_extension('ccomp', sources=['ccomp.c'],
include_dirs=[get_numpy_include_dirs()])
@@ -31,6 +32,8 @@ def configuration(parent_package='', top_path=None):
include_dirs=[get_numpy_include_dirs()])
config.add_extension('_convex_hull', sources=['_convex_hull.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('_greyreconstruct', sources=['_greyreconstruct.c'],
include_dirs=[get_numpy_include_dirs()])
return config
skimage/morphology/tests/test_reconstruction.py
+82
View File
@@ -0,0 +1,82 @@
"""
These tests are originally part of CellProfiler, code licensed under both GPL and BSD licenses.
Website: http://www.cellprofiler.org
Copyright (c) 2003-2009 Massachusetts Institute of Technology
Copyright (c) 2009-2011 Broad Institute
All rights reserved.
Original author: Lee Kamentsky
"""
import numpy as np
from numpy.testing import assert_array_almost_equal as assert_close
from skimage.morphology.greyreconstruct import reconstruction
def test_zeros():
"""Test reconstruction with image and mask of zeros"""
assert np.all(reconstruction(np.zeros((5, 7)), np.zeros((5, 7))) == 0)
def test_image_equals_mask():
"""Test reconstruction where the image and mask are the same"""
assert np.all(reconstruction(np.ones((7, 5)), np.ones((7, 5))) == 1)
def test_image_less_than_mask():
"""Test reconstruction where the image is uniform and less than mask"""
image = np.ones((5, 5))
mask = np.ones((5, 5)) * 2
assert np.all(reconstruction(image, mask) == 1)
def test_one_image_peak():
"""Test reconstruction with one peak pixel"""
image = np.ones((5, 5))
image[2, 2] = 2
mask = np.ones((5, 5)) * 3
assert np.all(reconstruction(image, mask) == 2)
def test_two_image_peaks():
"""Test reconstruction with two peak pixels isolated by the mask"""
image = np.array([[1, 1, 1, 1, 1, 1, 1, 1],
[1, 2, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 3, 1],
[1, 1, 1, 1, 1, 1, 1, 1]])
mask = np.array([[4, 4, 4, 1, 1, 1, 1, 1],
[4, 4, 4, 1, 1, 1, 1, 1],
[4, 4, 4, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 4, 4, 4],
[1, 1, 1, 1, 1, 4, 4, 4],
[1, 1, 1, 1, 1, 4, 4, 4]])
expected = np.array([[2, 2, 2, 1, 1, 1, 1, 1],
[2, 2, 2, 1, 1, 1, 1, 1],
[2, 2, 2, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 3, 3, 3],
[1, 1, 1, 1, 1, 3, 3, 3],
[1, 1, 1, 1, 1, 3, 3, 3]])
assert np.all(reconstruction(image, mask) == expected)
def test_zero_image_one_mask():
"""Test reconstruction with an image of all zeros and a mask that's not"""
result = reconstruction(np.zeros((10, 10)), np.ones((10, 10)))
assert np.all(result == 0)
def test_fill_hole():
"""Test reconstruction by erosion, which should fill holes in mask."""
seed = np.array([0, 8, 8, 8, 8, 8, 8, 8, 8, 0])
mask = np.array([0, 3, 6, 2, 1, 1, 1, 4, 2, 0])
result = reconstruction(seed, mask, method='erosion')
assert_close(result, np.array([0, 3, 6, 4, 4, 4, 4, 4, 2, 0]))
if __name__ == '__main__':
from numpy import testing
testing.run_module_suite()