diff --git a/skimage/morphology/_greyreconstruct.pyx b/skimage/morphology/_greyreconstruct.pyx index a59143f5..dff670bb 100644 --- a/skimage/morphology/_greyreconstruct.pyx +++ b/skimage/morphology/_greyreconstruct.pyx @@ -18,14 +18,14 @@ cimport cython @cython.boundscheck(False) def reconstruction_loop(np.ndarray[dtype=np.uint32_t, ndim=1, - negative_indices=False, mode='c'] avalues, + negative_indices=False, mode='c'] rank_array, 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, + np.int32_t current_idx, int image_stride): """The inner loop for reconstruction. @@ -41,66 +41,66 @@ def reconstruction_loop(np.ndarray[dtype=np.uint32_t, ndim=1, Parameters ---------- - avalues : array + rank_array : array The rank order of the flattened seed and mask images. aprev, anext: arrays Indices of previous and next pixels in rank sorted order. astrides : array Strides to neighbors of the current pixel. - current : int + current_idx : int Index of lowest-ranked pixel used as starting point in reconstruction loop. image_stride : int - Stride between seed image and mask image in `avalues`. + Stride between seed image and mask image in `rank_array`. """ cdef: - np.int32_t neighbor - np.uint32_t neighbor_value - np.uint32_t current_value - np.uint32_t mask_value + np.int32_t neighbor_idx + np.uint32_t neighbor_rank + np.uint32_t current_rank + np.uint32_t mask_rank np.int32_t current_link int i np.int32_t nprev np.int32_t nnext int nstrides = astrides.shape[0] - np.uint32_t *values = (avalues.data) + np.uint32_t *ranks = (rank_array.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: + 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 = current + strides[i] - neighbor_value = values[neighbor] + neighbor_idx = current_idx + strides[i] + neighbor_rank = ranks[neighbor_idx] # Only propagate neighbors ranked below the current rank - if neighbor_value < current_value: - mask_value = values[neighbor + image_stride] + if neighbor_rank < current_rank: + mask_rank = ranks[neighbor_idx + image_stride] # Only propagate neighbors ranked below the mask rank - if neighbor_value < mask_value: + if neighbor_rank < mask_rank: # Raise the neighbor to the mask rank if # the mask ranked below the current rank - if mask_value < current_value: - current_link = neighbor + image_stride - values[neighbor] = mask_value + if mask_rank < current_rank: + current_link = neighbor_idx + image_stride + ranks[neighbor_idx] = mask_rank else: - current_link = current - values[neighbor] = current_value + current_link = current_idx + ranks[neighbor_idx] = current_rank # unlink the neighbor - nprev = prev[neighbor] - nnext = next[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] = nnext - prev[neighbor] = current_link + next[neighbor_idx] = nnext + prev[neighbor_idx] = current_link if nnext >= 0: - prev[nnext] = neighbor - next[current_link] = neighbor - current = next[current] + 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 index d1434246..042d86f2 100644 --- a/skimage/morphology/greyreconstruct.py +++ b/skimage/morphology/greyreconstruct.py @@ -1,6 +1,6 @@ """ -`reconstruction` originally part of CellProfiler, code licensed under both GPL -and BSD licenses. +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 @@ -14,7 +14,7 @@ import numpy as np from skimage.filter.rank_order import rank_order -def reconstruction(image, mask, selem=None, offset=None, method='dilation'): +def reconstruction(seed, mask, selem=None, offset=None, method='dilation'): """Perform a morphological reconstruction of an image. Reconstruction requires a "seed" image and a "mask" image of equal shape. @@ -23,7 +23,7 @@ def reconstruction(image, mask, selem=None, offset=None, method='dilation'): Parameters ---------- - image : ndarray + seed : ndarray The seed image; a.k.a. marker image. mask : ndarray The maximum allowed value at each point. @@ -87,11 +87,11 @@ def reconstruction(image, mask, selem=None, offset=None, method='dilation'): transforms, but don't require a structuring element. """ - assert tuple(image.shape) == tuple(mask.shape) - if method == 'dilation' and np.any(image > mask): + 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(image < mask): + 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: @@ -100,7 +100,7 @@ def reconstruction(image, mask, selem=None, offset=None, method='dilation'): raise ImportError("_greyreconstruct extension not available.") if selem is None: - selem = np.ones([3]*image.ndim, dtype=bool) + selem = np.ones([3] * seed.ndim, dtype=bool) else: selem = selem.copy() @@ -113,54 +113,55 @@ def reconstruction(image, mask, selem=None, offset=None, method='dilation'): # 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 = 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(image) + pad_value = np.min(seed) elif method == 'erosion': - pad_value = np.max(image) - values = np.ones(dims) * pad_value - values[[0] + inside_slices] = image - values[[1] + inside_slices] = mask + 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(values.strides[1:]) / values.dtype.itemsize - image_stride = values.strides[0] / values.dtype.itemsize + 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) - values = values.flatten() - index_sorted = np.argsort(-values).astype(np.int32) - if method == 'erosion': + + 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(values), np.int32) - next = -np.ones(len(values), np.int32) + 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:] - # Create a rank-order value array so that the Cython inner-loop - # can operate on a uniform data type + # Cython inner-loop compares the rank of pixel values. if method == 'dilation': - value_rank, value_map = rank_order(values) + value_rank, value_map = rank_order(images) elif method == 'erosion': - value_rank, value_map = rank_order(-values) + value_rank, value_map = rank_order(-images) value_map = -value_map - current = index_sorted[0] - reconstruction_loop(value_rank, prev, next, nb_strides, current, - image_stride) + 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(image.shape) + 2 * padding + rec_img.shape = np.array(seed.shape) + 2 * padding return rec_img[inside_slices]