wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-04 17:23:16 +08:00
Code Issues Packages Projects Releases Wiki Activity

STY: Rename variables for clarity.

Browse Source 
In particular, it wasn't clear whether `image` was the seed image or the mask image rnd `values` was used to refer to both image intensity values and their rank-order.
This commit is contained in:
Tony S Yu
2012-08-18 19:21:55 -04:00
parent 7a56a7f35e
commit ab7626da3d
2 changed files with 62 additions and 61 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/morphology/_greyreconstruct.pyx
+31 -31
View File
@@ -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 = <np.uint32_t *>(avalues.data)
np.uint32_t *ranks = <np.uint32_t *>(rank_array.data)
np.int32_t *prev = <np.int32_t *>(aprev.data)
np.int32_t *next = <np.int32_t *>(anext.data)
np.int32_t *strides = <np.int32_t *>(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]
skimage/morphology/greyreconstruct.py
+31 -30
View File
@@ -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]