diff --git a/DEPENDS.txt b/DEPENDS.txt index b2858256..7793abdb 100644 --- a/DEPENDS.txt +++ b/DEPENDS.txt @@ -2,7 +2,7 @@ Build Requirements ------------------ * `Python >= 2.5 `__ * `Numpy >= 1.6 `__ -* `Cython >= 0.15 `__ +* `Cython >= 0.16 `__ `Matplotlib >= 1.0 `__ is needed to generate the examples in the documentation. diff --git a/doc/examples/plot_holes_and_peaks.py b/doc/examples/plot_holes_and_peaks.py new file mode 100644 index 00000000..a9c3a7fe --- /dev/null +++ b/doc/examples/plot_holes_and_peaks.py @@ -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 +""" diff --git a/doc/examples/plot_regional_maxima.py b/doc/examples/plot_regional_maxima.py new file mode 100644 index 00000000..9d4de9b1 --- /dev/null +++ b/doc/examples/plot_regional_maxima.py @@ -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. +""" diff --git a/skimage/morphology/__init__.py b/skimage/morphology/__init__.py index abc9986f..efa2077d 100644 --- a/skimage/morphology/__init__.py +++ b/skimage/morphology/__init__.py @@ -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 diff --git a/skimage/morphology/_greyreconstruct.pyx b/skimage/morphology/_greyreconstruct.pyx new file mode 100644 index 00000000..5d2cb0c5 --- /dev/null +++ b/skimage/morphology/_greyreconstruct.pyx @@ -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] + diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py new file mode 100644 index 00000000..9362e3fb --- /dev/null +++ b/skimage/morphology/greyreconstruct.py @@ -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] + diff --git a/skimage/morphology/setup.py b/skimage/morphology/setup.py index 2dd4a378..f5dd7f19 100644 --- a/skimage/morphology/setup.py +++ b/skimage/morphology/setup.py @@ -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 diff --git a/skimage/morphology/tests/test_reconstruction.py b/skimage/morphology/tests/test_reconstruction.py new file mode 100644 index 00000000..74f61b9f --- /dev/null +++ b/skimage/morphology/tests/test_reconstruction.py @@ -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()