diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py index 04ea54a3..9362e3fb 100644 --- a/skimage/morphology/greyreconstruct.py +++ b/skimage/morphology/greyreconstruct.py @@ -14,26 +14,41 @@ import numpy as np from skimage.filter.rank_order import rank_order -def reconstruction(seed, mask, selem=None, offset=None, method='dilation'): +def reconstruction(seed, mask, method='dilation', selem=None, offset=None): """Perform a morphological reconstruction of an image. - Reconstruction requires a "seed" image and a "mask" image of equal shape. - These images set the minimum and maximum possible values of the - reconstructed 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. + The seed image (a.k.a. marker image), which specifies the values that + are dilated or eroded. mask : ndarray - The maximum allowed value at each point. + 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. - method : {'dilation'|'erosion'} - Perform reconstruction by dilation or erosion. In dilation (erosion), - the seed image is dilated (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. Returns ------- @@ -42,11 +57,29 @@ def reconstruction(seed, mask, selem=None, offset=None, method='dilation'): Examples -------- - Here, we try to extract the bright features of an image by subtracting a - background image created by reconstruction. - >>> 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)