From 33a19a8c7bb8560af53cc7b0d91e675862c6a8ca Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Fri, 17 Aug 2012 22:48:02 -0400 Subject: [PATCH] DOC: add comments to clarify algorithm --- skimage/morphology/greyreconstruct.py | 39 ++++++++++++++++----------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py index 145d007a..8c12a9c5 100644 --- a/skimage/morphology/greyreconstruct.py +++ b/skimage/morphology/greyreconstruct.py @@ -1,5 +1,6 @@ """ -`reconstruction` originally part of CellProfiler, code licensed under both GPL and BSD licenses. +`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 @@ -14,17 +15,18 @@ from skimage.filter.rank_order import rank_order def reconstruction(image, mask, selem=None, offset=None): - """Perform a morphological reconstruction of the image. + """Perform a morphological reconstruction of an 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" + Reconstruction requires a "seed" image and a "mask" image. Currently, this + only implements reconstruction by dilation, such that the seed image is + dilated until it is constrained by the mask. Thus, he "seed" and "mask" images will be the minimum and maximum possible values of the reconstructed image, respectively. Parameters ---------- image : ndarray - The seed image. + The seed image; a.k.a. marker image. mask : ndarray The maximum allowed value at each point. selem : ndarray @@ -42,9 +44,13 @@ def reconstruction(image, mask, selem=None, offset=None): 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) + + [1] Vincent, L., "Morphological Grayscale Reconstruction in Image Analysis: + Applications and Efficient Algorithms", IEEE Transactions on Image + Processing (1993) + + [2] Soille, P., "Morphological Image Analysis: Principles and Applications", + Chapter 6, 2nd edition (2003), ISBN 3540429883. Examples -------- @@ -98,27 +104,27 @@ def reconstruction(image, mask, selem=None, offset=None): # 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 + # Make padding for edges of reconstructed image so we can ignore boundaries padding = (np.array(selem.shape) / 2).astype(int) dims = np.zeros(image.ndim + 1, dtype=int) dims[1:] = np.array(image.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. 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 + # 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) + nb_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) @@ -130,10 +136,11 @@ def reconstruction(image, mask, selem=None, offset=None): # Create a rank-order value array so that the Cython inner-loop # can operate on a uniform data type + # fragile: `reconstruction_loop` needs 'uint32' conversion by `rank_order` values, value_map = rank_order(values) current = value_sort[0] - reconstruction_loop(values, prev, next, strides, current, image_stride) + reconstruction_loop(values, prev, next, nb_strides, current, image_stride) # Reshape the values array to the shape of the padded image # and return the unpadded portion of that result