diff --git a/doc/examples/applications/plot_peak_detection_comparison.py b/doc/examples/applications/plot_peak_detection_comparison.py new file mode 100644 index 00000000..2cf7f933 --- /dev/null +++ b/doc/examples/applications/plot_peak_detection_comparison.py @@ -0,0 +1,218 @@ +""" +============== +Peak detection +============== + +Peak detection (a.k.a. spot detection or particle detection) is a common +image analysis step. For example, it's used to detect tracer particles in a +flow for particle image velocimetry (i.e. PIV) and to identify features in +the Hough transform. + +To simplify plotting code, let's define a simple function that creates a new +figure on each call, and removes tick labels. +""" + +import matplotlib.pyplot as plt + +plt.rcParams['axes.titlesize'] = 10 +plt.rcParams['font.size'] = 10 +def imshow(image, **kwargs): + plt.figure(figsize=(2.5, 2.5)) + plt.imshow(image, **kwargs) + plt.axis('off') + +""" +To explore different peak detection techniques, we use an image of circles with +added noise: +""" + +from skimage import data +img = data.load('noisy_circles.jpg') +imshow(img) + +""" + +.. image:: PLOT2RST.current_figure + +This image is noisy and has uneven background illumination. The peaks in the +image, while readily identified by eye, can be tricky to find algorithmically. +The first thing we need to do is remove the high-frequency noise; this can +be done with a simple Gaussian filter. +""" + +import scipy.ndimage as ndimg +img_smooth = ndimg.gaussian_filter(img, 3) + +imshow(img_smooth) + +""" + +.. image:: PLOT2RST.current_figure + +Thresholding +============ + +One way to extract the background is to threshold the image. +""" + +thresh_value = 100 +background = img_smooth.copy() +background[img_smooth > thresh_value] = 0 +peaks = img_smooth - background + +""" +Here, all pixels values below the threshold value are subtracted from the +image. The resulting background image and the extracted peaks are shown below. +""" + +imshow(background, vmin=0, vmax=255) +plt.title("background image (thresholding)") + +""" +.. image:: PLOT2RST.current_figure +""" + +imshow(peaks, vmin=0, vmax=255) +plt.title("peaks (thresholding)") + +""" +.. image:: PLOT2RST.current_figure + +Because of uneven illumination, peaks on the right bleed into each other. +Increasing the threshold will fix this problem, but it will also cause some +peaks on the left to go undetected. + +Morphological reconstruction +============================ +""" + +import numpy as np +img_r = np.int32(img_smooth) + +import skimage.morphology as morph +h = 20 +rec = morph.reconstruction(img_r-h, img_r) + +imshow(img_r, vmin=0, vmax=255) +plt.title("original (smoothed) image") + +""" +.. image:: PLOT2RST.current_figure +""" + +imshow(rec, vmin=0, vmax=255) +plt.title("background image (reconstruction)") + +""" +.. image:: PLOT2RST.current_figure + +This reconstructed image looks pretty much like the original, except that the +peaks in the image are truncated. The reconstructed image can then be +subtracted from the original image to reveal the peaks of the image. +""" + +imshow(img_r-rec) +plt.title("h-dome of image") + +""" +.. image:: PLOT2RST.current_figure + +The result is known as the h-dome transformation [2]_, which extracts peaks of +height `h` from the original image. To better understand what's going on, +let's take a 1D slice along the middle of the image (cutting through peaks in +the image). +""" + +img_slice = img_r[99:100, :] +rec_slice = morph.reconstruction(img_slice-h, img_slice) + +""" +Plotting the reconstructed image (slice) next to the original image and the +seed image shed light on the reconstruction process +""" +plt.figure(figsize=(4, 3)) +plt.plot(img_slice[0], 'k', label='original image') +plt.plot(img_slice[0]-h, '0.5', label='seed image') +plt.plot(rec_slice[0], 'r', label='reconstructed') +plt.title("image slice") +plt.xlabel('x') +plt.ylabel('intensity') +plt.legend() + +""" +.. image:: PLOT2RST.current_figure + +Here, you see that morphological reconstruction dilates the seed image (i.e. +the `h`-shifted image) until it intersects the mask (original image). Note that +the peaks in the original image have very different intensity values (e.g. the +peak at x=200 and x=100 differ by about 80). Subtracting the reconstructed +image from the original image gives peaks of roughly equal intensity. Thus, the +h-dome transformation is quiet effective at removing uneven, dark backgrounds +from bright features. The inverse operation---the h-basin +transformation---should be used when removing bright backgrounds from dark +features. + + +White tophat +============ +""" +selem = morph.disk(10) +img_t = np.uint8(img_smooth) +opening = morph.greyscale_open(img_t, selem) +top_hat = img_t - opening + +imshow(opening, vmin=0, vmax=255) +plt.title("Greyscale opening of image") + +""" +.. image:: PLOT2RST.current_figure +""" + + +imshow(top_hat) +plt.title("Tophat with disk of r = 10") + +""" +.. image:: PLOT2RST.current_figure +""" + +selem = morph.disk(5) +top_hat = morph.greyscale_white_top_hat(img_t, selem) + +imshow(top_hat) +plt.title("Tophat with disk of r = 5") + +""" +.. image:: PLOT2RST.current_figure +""" + +selem = morph.square(20) +opening = morph.greyscale_open(img_t, selem) +# scikit's top hat filter uses uint8 and doesn't check for over(under)flow. +mask = opening > img_t +opening[mask] = img_t[mask] +top_hat = img_t - opening + +imshow(opening, vmin=0, vmax=255) +plt.title("Greyscale opening of image") + +""" +.. image:: PLOT2RST.current_figure +""" + +imshow(top_hat) +plt.title("Tophat with square of w = 10") + +plt.show() + +""" +.. image:: PLOT2RST.current_figure + + +References +========== + +.. [1] Crocker and Grier, Journal of Colloid and Interface Science (1996) +.. [2] Vincent, L., IEEE Transactions on Image Processing (1993) + +""" diff --git a/skimage/data/noisy_circles.jpg b/skimage/data/noisy_circles.jpg new file mode 100644 index 00000000..a0f49d67 Binary files /dev/null and b/skimage/data/noisy_circles.jpg differ diff --git a/skimage/morphology/__init__.py b/skimage/morphology/__init__.py index d639be1c..24982269 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..18197812 --- /dev/null +++ b/skimage/morphology/_greyreconstruct.pyx @@ -0,0 +1,85 @@ +""" +`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 + +""" + +from __future__ import division +import numpy as np + +cimport numpy as np +cimport cython + + +@cython.boundscheck(False) +def reconstruction_loop(np.ndarray[dtype=np.uint32_t, ndim=1, + negative_indices = False, + mode = 'c'] avalues, + np.ndarray[dtype=np.int32_t, ndim=1, + negative_indices = False, + mode = 'c'] aprev, + np.ndarray[dtype=np.int32_t, ndim=1, + negative_indices = False, + mode = 'c'] anext, + np.ndarray[dtype=np.int32_t, ndim=1, + negative_indices = False, + mode = 'c'] astrides, + np.int32_t current, + int image_stride): + """The inner loop for reconstruction""" + cdef: + np.int32_t neighbor + np.uint32_t neighbor_value + np.uint32_t current_value + np.uint32_t mask_value + np.int32_t link + int i + np.int32_t nprev + np.int32_t nnext + int nstrides = astrides.shape[0] + np.uint32_t *values = (avalues.data) + np.int32_t *prev = (aprev.data) + np.int32_t *next = (anext.data) + np.int32_t *strides = (astrides.data) + + while current != -1: + if current < image_stride: + current_value = values[current] + if current_value == 0: + break + for i in range(nstrides): + neighbor = current + strides[i] + neighbor_value = values[neighbor] + # Only do neighbors less than the current value + if neighbor_value < current_value: + mask_value = values[neighbor + image_stride] + # Only do neighbors less than the mask value + if neighbor_value < mask_value: + # Raise the neighbor to the mask value if + # the mask is less than current + if mask_value < current_value: + link = neighbor + image_stride + values[neighbor] = mask_value + else: + link = current + values[neighbor] = current_value + # unlink the neighbor + nprev = prev[neighbor] + nnext = next[neighbor] + next[nprev] = nnext + if nnext != -1: + prev[nnext] = nprev + # link the neighbor after the link + nnext = next[link] + next[neighbor] = nnext + prev[neighbor] = link + if nnext >= 0: + prev[nnext] = neighbor + next[link] = neighbor + current = next[current] + diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py new file mode 100644 index 00000000..7f0def38 --- /dev/null +++ b/skimage/morphology/greyreconstruct.py @@ -0,0 +1,150 @@ +""" +`reconstruction` 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 skimage.filter.rank_order import rank_order + + +def reconstruction(image, mask, selem=None, offset=None): + """Perform a morphological reconstruction of the image. + + Reconstruction requires a "seed" image and a "mask" image. The seed image + gets dilated until it is constrained by the mask. The "seed" and "mask" + images will be the minimum and maximum possible values of the reconstructed + image. + + Parameters + ---------- + image : ndarray + The seed image. + + mask : ndarray + The maximum allowed value at each point. + + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + + Returns + ------- + reconstructed : ndarray + The result of morphological reconstruction. + + Notes + ----- + The algorithm is taken from: + Robinson, "Efficient morphological reconstruction: a downhill filter", + Pattern Recognition Letters 25 (2004) 1759-1767. + + Applications for greyscale reconstruction are discussed in: + Vincent, L., "Morphological Grayscale Reconstruction in Image Analysis: + Applications and Efficient Algorithms", IEEE Transactions on Image + Processing (1993) + + Examples + -------- + Uses for greyscale reconstruction are described in Vincent (1993). For + example, let's try to extract the features of an image by subtracting a + background image created by reconstruction. + + First, create an image where the "bumps" are the features that + we want to extract: + + >>> import numpy as np + >>> from scikits.image.morphology.grey import grey_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 + >>> rec = grey_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 - rec + + This operation is known as the h-dome of the image, which leaves features + of height `h` in the subtracted image. The h-dome transform, and its + inverse h-basin, are analogous to the white top-hat and black top-hat + transforms, but don't require a structuring element. + + """ + assert tuple(image.shape) == tuple(mask.shape) + assert np.all(image <= mask) + try: + from ._morphrec import reconstruction_loop + except ImportError: + raise ImportError("_morphrec extension not available.") + + if selem is None: + selem = np.ones([3]*image.ndim, bool) + else: + selem = selem.copy() + + if offset == None: + assert all([d % 2 == 1 for d in selem.shape]),\ + "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 + # + # Construct an array that's padded on the edges so we can ignore boundaries + # The array is a dstack of the image and the mask; this lets us interleave + # image and mask pixels when sorting which makes list manipulations easier + # + padding = (np.array(selem.shape)/2).astype(int) + dims = np.zeros(image.ndim+1,int) + dims[1:] = np.array(image.shape)+2*padding + dims[0] = 2 + inside_slices = [slice(p,-p) for p in padding] + values = np.ones(dims)*np.min(image) + values[[0]+inside_slices] = image + values[[1]+inside_slices] = mask + # + # Create a list of strides across the array to get the neighbors + # within a flattened array + # + value_stride = np.array(values.strides[1:]) / values.dtype.itemsize + image_stride = values.strides[0] / values.dtype.itemsize + selem_mgrid = np.mgrid[[slice(-o,d - o) + for d,o in zip(selem.shape,offset)]] + selem_offsets = selem_mgrid[:,selem].transpose() + strides = np.array([np.sum(value_stride * selem_offset) + for selem_offset in selem_offsets], np.int32) + values = values.flatten() + value_sort = np.lexsort([-values]).astype(np.int32) + # + # Make a linked list of pixels sorted by value. -1 is the list terminator. + # + prev = -np.ones(len(values), np.int32) + next = -np.ones(len(values), np.int32) + prev[value_sort[1:]] = value_sort[:-1] + next[value_sort[:-1]] = value_sort[1:] + # + # Create a rank-order value array so that the Cython inner-loop + # can operate on a uniform data type + # + values, value_map = rank_order(values) + current = value_sort[0] + + reconstruction_loop(values, prev, next, strides, current, image_stride) + # + # Reshape the values array to the shape of the padded image + # and return the unpadded portion of that result + # + values = value_map[values[:image_stride]] + values.shape = np.array(image.shape)+2*padding + return values[inside_slices] + diff --git a/skimage/morphology/setup.py b/skimage/morphology/setup.py index fcf33f7f..fda13679 100644 --- a/skimage/morphology/setup.py +++ b/skimage/morphology/setup.py @@ -18,6 +18,7 @@ def configuration(parent_package='', top_path=None): cython(['_skeletonize.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..61b72033 --- /dev/null +++ b/skimage/morphology/tests/test_reconstruction.py @@ -0,0 +1,73 @@ +""" +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 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) + + +if __name__ == '__main__': + from numpy import testing + testing.run_module_suite()