mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-09 12:51:42 +08:00
Merge Emmanuelle's port of CellProfiler's medial axis skeletonization.
This commit is contained in:
Stefan van der Walt
@@ -0,0 +1,70 @@
|
||||
"""
|
||||
===========================
|
||||
Medial axis skeletonization
|
||||
===========================
|
||||
|
||||
The medial axis of an object is the set of all points having more than one
|
||||
closest point on the object's boundary. It is often called the **topological
|
||||
skeleton**, because it is a 1-pixel wide skeleton of the object, with the same
|
||||
connectivity as the original object.
|
||||
|
||||
Here, we use the medial axis transform to compute the width of the foreground
|
||||
objects. As the function ``medial_axis`` (``skimage.morphology.medial_axis``)
|
||||
returns the distance transform in addition to the medial axis (with the keyword
|
||||
argument ``return_distance=True``), it is possible to compute the distance to
|
||||
the background for all points of the medial axis with this function. This gives
|
||||
an estimate of the local width of the objects.
|
||||
|
||||
For a skeleton with fewer branches, there exists another skeletonization
|
||||
algorithm in ``skimage``: ``skimage.morphology.skeletonize``, that computes
|
||||
a skeleton by iterative morphological thinnings.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
from scipy import ndimage
|
||||
from skimage.morphology import medial_axis
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
|
||||
def microstructure(l=256):
|
||||
"""
|
||||
Synthetic binary data: binary microstructure with blobs.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
|
||||
l: int, optional
|
||||
linear size of the returned image
|
||||
|
||||
"""
|
||||
n = 5
|
||||
x, y = np.ogrid[0:l, 0:l]
|
||||
mask_outer = (x - l/2)**2 + (y - l/2)**2 < (l/2)**2
|
||||
mask = np.zeros((l, l))
|
||||
generator = np.random.RandomState(1)
|
||||
points = l * generator.rand(2, n**2)
|
||||
mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
|
||||
mask = ndimage.gaussian_filter(mask, sigma=l/(4.*n))
|
||||
return mask > mask.mean()
|
||||
|
||||
data = microstructure(l=64)
|
||||
|
||||
# Compute the medial axis (skeleton) and the distance transform
|
||||
skel, distance = medial_axis(data, return_distance=True)
|
||||
|
||||
# Distance to the background for pixels of the skeleton
|
||||
dist_on_skel = distance * skel
|
||||
|
||||
plt.figure(figsize=(8, 4))
|
||||
plt.subplot(121)
|
||||
plt.imshow(data, cmap=plt.cm.gray, interpolation='nearest')
|
||||
plt.axis('off')
|
||||
plt.subplot(122)
|
||||
plt.imshow(dist_on_skel, cmap=plt.cm.spectral, interpolation='nearest')
|
||||
plt.contour(data, [0.5], colors='w')
|
||||
plt.axis('off')
|
||||
|
||||
plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
|
||||
right=1)
|
||||
plt.show()
|
||||
|
||||
@@ -15,8 +15,8 @@ results. The input is a 2D ndarray, with either boolean or integer elements.
|
||||
In the case of boolean, 'True' indicates foreground, and for integer arrays,
|
||||
the foreground is 1's.
|
||||
"""
|
||||
from scikits.image.morphology import skeletonize
|
||||
from scikits.image.draw import draw
|
||||
from skimage.morphology import skeletonize
|
||||
from skimage.draw import draw
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
|
||||
@@ -2,4 +2,4 @@ from grey import *
|
||||
from selem import *
|
||||
from .ccomp import label
|
||||
from watershed import watershed, is_local_maximum
|
||||
from skeletonize import skeletonize
|
||||
from skeletonize import skeletonize, medial_axis
|
||||
|
||||
@@ -0,0 +1,212 @@
|
||||
'''
|
||||
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
|
||||
cimport numpy as np
|
||||
cimport cython
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
def _skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2,
|
||||
negative_indices=False, mode='c'] result,
|
||||
np.ndarray[dtype=np.int32_t, ndim=1,
|
||||
negative_indices=False, mode='c'] i,
|
||||
np.ndarray[dtype=np.int32_t, ndim=1,
|
||||
negative_indices=False, mode='c'] j,
|
||||
np.ndarray[dtype=np.int32_t, ndim=1,
|
||||
negative_indices=False, mode='c'] order,
|
||||
np.ndarray[dtype=np.uint8_t, ndim=1,
|
||||
negative_indices=False, mode='c'] table):
|
||||
"""
|
||||
Inner loop of skeletonize function
|
||||
|
||||
Parameters
|
||||
----------
|
||||
|
||||
result : ndarray of uint8
|
||||
On input, the image to be skeletonized, on output the skeletonized
|
||||
image.
|
||||
|
||||
i, j : ndarrays
|
||||
The coordinates of each foreground pixel in the image
|
||||
|
||||
order : ndarray
|
||||
The index of each pixel, in the order of processing (order[0] is
|
||||
the first pixel to process, etc.)
|
||||
|
||||
table : ndarray
|
||||
The 512-element lookup table of values after transformation
|
||||
(whether to keep or not each configuration in a binary 3x3 array)
|
||||
|
||||
Notes
|
||||
-----
|
||||
|
||||
The loop determines whether each pixel in the image can be removed without
|
||||
changing the Euler number of the image. The pixels are ordered by
|
||||
increasing distance from the background which means a point nearer to
|
||||
the quench-line of the brushfire will be evaluated later than a
|
||||
point closer to the edge.
|
||||
|
||||
Note that the neighbourhood of a pixel may evolve before the loop
|
||||
arrives at this pixel. This is why it is possible to compute the
|
||||
skeleton in only one pass, thanks to an adapted ordering of the
|
||||
pixels.
|
||||
"""
|
||||
cdef:
|
||||
np.int32_t accumulator
|
||||
np.int32_t index, order_index
|
||||
np.int32_t ii, jj
|
||||
|
||||
for index in range(order.shape[0]):
|
||||
accumulator = 16
|
||||
order_index = order[index]
|
||||
ii = i[order_index]
|
||||
jj = j[order_index]
|
||||
# Compute the configuration around the pixel
|
||||
if ii > 0:
|
||||
if jj > 0 and result[ii - 1, jj - 1]:
|
||||
accumulator += 1
|
||||
if result[ii - 1, jj]:
|
||||
accumulator += 2
|
||||
if jj < result.shape[1] - 1 and result[ii - 1, jj + 1]:
|
||||
accumulator += 4
|
||||
if jj > 0 and result[ii, jj - 1]:
|
||||
accumulator += 8
|
||||
if jj < result.shape[1] - 1 and result[ii, jj + 1]:
|
||||
accumulator += 32
|
||||
if ii < result.shape[0]-1:
|
||||
if jj > 0 and result[ii + 1, jj - 1]:
|
||||
accumulator += 64
|
||||
if result[ii + 1, jj]:
|
||||
accumulator += 128
|
||||
if jj < result.shape[1] - 1 and result[ii + 1, jj + 1]:
|
||||
accumulator += 256
|
||||
# Assign the value of table corresponding to the configuration
|
||||
result[ii, jj] = table[accumulator]
|
||||
|
||||
@cython.boundscheck(False)
|
||||
def _table_lookup_index(np.ndarray[dtype=np.uint8_t, ndim=2,
|
||||
negative_indices=False, mode='c'] image):
|
||||
"""
|
||||
Return an index into a table per pixel of a binary image
|
||||
|
||||
Take the sum of true neighborhood pixel values where the neighborhood
|
||||
looks like this::
|
||||
|
||||
1 2 4
|
||||
8 16 32
|
||||
64 128 256
|
||||
|
||||
This code could be replaced by a convolution with the kernel::
|
||||
|
||||
256 128 64
|
||||
32 16 8
|
||||
4 2 1
|
||||
|
||||
but this runs about twice as fast because of inlining and the
|
||||
hardwired kernel.
|
||||
"""
|
||||
cdef:
|
||||
np.ndarray[dtype=np.int32_t, ndim=2,
|
||||
negative_indices=False, mode='c'] indexer
|
||||
np.int32_t *p_indexer
|
||||
np.uint8_t *p_image
|
||||
np.int32_t i_stride
|
||||
np.int32_t i_shape
|
||||
np.int32_t j_shape
|
||||
np.int32_t i
|
||||
np.int32_t j
|
||||
np.int32_t offset
|
||||
|
||||
i_shape = image.shape[0]
|
||||
j_shape = image.shape[1]
|
||||
indexer = np.zeros((i_shape, j_shape), np.int32)
|
||||
p_indexer = <np.int32_t *>indexer.data
|
||||
p_image = <np.uint8_t *>image.data
|
||||
i_stride = image.strides[0]
|
||||
assert i_shape >= 3 and j_shape >= 3, \
|
||||
"Please use the slow method for arrays < 3x3"
|
||||
|
||||
for i in range(1, i_shape-1):
|
||||
offset = i_stride* i + 1
|
||||
for j in range(1, j_shape - 1):
|
||||
if p_image[offset]:
|
||||
p_indexer[offset + i_stride + 1] += 1
|
||||
p_indexer[offset + i_stride] += 2
|
||||
p_indexer[offset + i_stride - 1] += 4
|
||||
p_indexer[offset + 1] += 8
|
||||
p_indexer[offset] += 16
|
||||
p_indexer[offset - 1] += 32
|
||||
p_indexer[offset - i_stride + 1] += 64
|
||||
p_indexer[offset - i_stride] += 128
|
||||
p_indexer[offset - i_stride - 1] += 256
|
||||
offset += 1
|
||||
#
|
||||
# Do the corner cases (literally)
|
||||
#
|
||||
if image[0, 0]:
|
||||
indexer[0, 0] += 16
|
||||
indexer[0, 1] += 8
|
||||
indexer[1, 0] += 2
|
||||
indexer[1, 1] += 1
|
||||
|
||||
if image[0, j_shape - 1]:
|
||||
indexer[0, j_shape - 2] += 32
|
||||
indexer[0, j_shape - 1] += 16
|
||||
indexer[1, j_shape - 2] += 4
|
||||
indexer[1, j_shape - 1] += 2
|
||||
|
||||
if image[i_shape - 1, 0]:
|
||||
indexer[i_shape - 2, 0] += 128
|
||||
indexer[i_shape - 2, 1] += 64
|
||||
indexer[i_shape - 1, 0] += 16
|
||||
indexer[i_shape - 1, 1] += 8
|
||||
|
||||
if image[i_shape - 1, j_shape - 1]:
|
||||
indexer[i_shape - 2, j_shape - 2] += 256
|
||||
indexer[i_shape - 2, j_shape - 1] += 128
|
||||
indexer[i_shape - 1, j_shape - 2] += 32
|
||||
indexer[i_shape - 1, j_shape - 1] += 16
|
||||
#
|
||||
# Do the edges
|
||||
#
|
||||
for j in range(1, j_shape - 1):
|
||||
if image[0, j]:
|
||||
indexer[0, j - 1] += 32
|
||||
indexer[0, j] += 16
|
||||
indexer[0, j + 1] += 8
|
||||
indexer[1, j - 1] += 4
|
||||
indexer[1, j] += 2
|
||||
indexer[1, j + 1] += 1
|
||||
if image[i_shape - 1, j]:
|
||||
indexer[i_shape - 2, j - 1] += 256
|
||||
indexer[i_shape - 2, j] += 128
|
||||
indexer[i_shape - 2, j + 1] += 64
|
||||
indexer[i_shape - 1, j - 1] += 32
|
||||
indexer[i_shape - 1, j] += 16
|
||||
indexer[i_shape - 1, j + 1] += 8
|
||||
|
||||
for i in range(1, i_shape - 1):
|
||||
if image[i, 0]:
|
||||
indexer[i - 1, 0] += 128
|
||||
indexer[i, 0] += 16
|
||||
indexer[i + 1, 0] += 2
|
||||
indexer[i - 1, 1] += 64
|
||||
indexer[i, 1] += 8
|
||||
indexer[i + 1, 1] += 1
|
||||
if image[i, j_shape - 1]:
|
||||
indexer[i - 1, j_shape - 2] += 256
|
||||
indexer[i, j_shape - 2] += 32
|
||||
indexer[i + 1, j_shape - 2] += 4
|
||||
indexer[i - 1, j_shape - 1] += 128
|
||||
indexer[i, j_shape - 1] += 16
|
||||
indexer[i + 1, j_shape - 1] += 2
|
||||
return indexer
|
||||
@@ -14,6 +14,7 @@ def configuration(parent_package='', top_path=None):
|
||||
cython(['ccomp.pyx'], working_path=base_path)
|
||||
cython(['cmorph.pyx'], working_path=base_path)
|
||||
cython(['_watershed.pyx'], working_path=base_path)
|
||||
cython(['_skeletonize.pyx'], working_path=base_path)
|
||||
|
||||
config.add_extension('ccomp', sources=['ccomp.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
@@ -21,6 +22,9 @@ def configuration(parent_package='', top_path=None):
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_watershed', sources=['_watershed.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_skeletonize', sources=['_skeletonize.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
|
||||
return config
|
||||
|
||||
|
||||
@@ -1,10 +1,13 @@
|
||||
"""Use an iterative thinning algorithm to find the skeletons of binary
|
||||
objects in an image.
|
||||
|
||||
"""
|
||||
Algorithms for computing the skeleton of a binary image
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
from scipy.ndimage import correlate
|
||||
from scipy import ndimage
|
||||
|
||||
from _skeletonize import _skeletonize_loop, _table_lookup_index
|
||||
|
||||
# --------- Skeletonization by morphological thinning ---------
|
||||
|
||||
def skeletonize(image):
|
||||
"""Return the skeleton of a binary image.
|
||||
@@ -24,6 +27,10 @@ def skeletonize(image):
|
||||
skeleton : ndarray
|
||||
A matrix containing the thinned image.
|
||||
|
||||
See also
|
||||
--------
|
||||
medial_axis
|
||||
|
||||
Notes
|
||||
-----
|
||||
The algorithm [1] works by making successive passes of the image,
|
||||
@@ -107,7 +114,7 @@ def skeletonize(image):
|
||||
pixelRemoved = False;
|
||||
|
||||
# assign each pixel a unique value based on its foreground neighbours
|
||||
neighbours = correlate(skeleton, mask, mode='constant')
|
||||
neighbours = ndimage.correlate(skeleton, mask, mode='constant')
|
||||
|
||||
# ignore background
|
||||
neighbours *= skeleton
|
||||
@@ -126,7 +133,7 @@ def skeletonize(image):
|
||||
skeleton[code_mask] = 0
|
||||
|
||||
# pass 2 - remove the 2's and 3's
|
||||
neighbours = correlate(skeleton, mask, mode='constant')
|
||||
neighbours = ndimage.correlate(skeleton, mask, mode='constant')
|
||||
neighbours *= skeleton
|
||||
codes = np.take(lut, neighbours)
|
||||
code_mask = (codes == 2)
|
||||
@@ -139,3 +146,226 @@ def skeletonize(image):
|
||||
skeleton[code_mask] = 0
|
||||
|
||||
return skeleton
|
||||
|
||||
# --------- Skeletonization by medial axis transform --------
|
||||
|
||||
_eight_connect = ndimage.generate_binary_structure(2, 2)
|
||||
|
||||
|
||||
def medial_axis(image, mask=None, return_distance=False):
|
||||
"""
|
||||
Compute the medial axis transform of a binary image
|
||||
|
||||
Parameters
|
||||
----------
|
||||
|
||||
image : binary ndarray
|
||||
|
||||
mask : binary ndarray, optional
|
||||
If a mask is given, only those elements with a true value in `mask`
|
||||
are used for computing the medial axis.
|
||||
|
||||
return_distance : bool, optional
|
||||
If true, the distance transform is returned as well as the skeleton.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
||||
out : ndarray of bools
|
||||
Medial axis transform of the image
|
||||
|
||||
dist : ndarray of ints
|
||||
Distance transform of the image (only returned if `return_distance`
|
||||
is True)
|
||||
|
||||
See also
|
||||
--------
|
||||
skeletonize
|
||||
|
||||
Notes
|
||||
-----
|
||||
This algorithm computes the medial axis transform of an image
|
||||
as the ridges of its distance transform.
|
||||
|
||||
The different steps of the algorithm are as follows
|
||||
* A lookup table is used, that assigns 0 or 1 to each configuration of
|
||||
the 3x3 binary square, whether the central pixel should be removed
|
||||
or kept. We want a point to be removed if it has more than one neighbor
|
||||
and if removing it does not change the number of connected components.
|
||||
|
||||
* The distance transform to the background is computed, as well as
|
||||
the cornerness of the pixel.
|
||||
|
||||
* The foreground (value of 1) points are ordered by
|
||||
the distance transform, then the cornerness.
|
||||
|
||||
* A cython function is called to reduce the image to its skeleton. It
|
||||
processes pixels in the order determined at the previous step, and
|
||||
removes or maintains a pixel according to the lookup table. Because
|
||||
of the ordering, it is possible to process all pixels in only one
|
||||
pass.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> square = np.zeros((7, 7), dtype=np.uint8)
|
||||
>>> square[1:-1, 2:-2] = 1
|
||||
>>> square
|
||||
array([[0, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 1, 1, 1, 0, 0],
|
||||
[0, 0, 1, 1, 1, 0, 0],
|
||||
[0, 0, 1, 1, 1, 0, 0],
|
||||
[0, 0, 1, 1, 1, 0, 0],
|
||||
[0, 0, 1, 1, 1, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
|
||||
>>> morphology.medial_axis(square).astype(np.uint8)
|
||||
array([[0, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 1, 0, 1, 0, 0],
|
||||
[0, 0, 0, 1, 0, 0, 0],
|
||||
[0, 0, 0, 1, 0, 0, 0],
|
||||
[0, 0, 0, 1, 0, 0, 0],
|
||||
[0, 0, 1, 0, 1, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
|
||||
|
||||
"""
|
||||
global _eight_connect
|
||||
if mask is None:
|
||||
masked_image = image.astype(np.bool)
|
||||
else:
|
||||
masked_image = image.astype(bool).copy()
|
||||
masked_image[~mask] = False
|
||||
#
|
||||
# Build lookup table - three conditions
|
||||
# 1. Keep only positive pixels (center_is_foreground array).
|
||||
# AND
|
||||
# 2. Keep if removing the pixel results in a different connectivity
|
||||
# (if the number of connected components is different with and
|
||||
# without the central pixel)
|
||||
# OR
|
||||
# 3. Keep if # pixels in neighbourhood is 2 or less
|
||||
# Note that table is independent of image
|
||||
center_is_foreground = (np.arange(512) & 2**4).astype(bool)
|
||||
table = (center_is_foreground # condition 1.
|
||||
&
|
||||
(np.array([ndimage.label(_pattern_of(index), _eight_connect)[1] !=
|
||||
ndimage.label(_pattern_of(index & ~ 2**4),
|
||||
_eight_connect)[1]
|
||||
for index in range(512)]) # condition 2
|
||||
|
|
||||
np.array([np.sum(_pattern_of(index)) < 3 for index in range(512)]))
|
||||
# condition 3
|
||||
)
|
||||
|
||||
|
||||
# Build distance transform
|
||||
distance = ndimage.distance_transform_edt(masked_image)
|
||||
if return_distance:
|
||||
store_distance = distance.copy()
|
||||
|
||||
# Corners
|
||||
# The processing order along the edge is critical to the shape of the
|
||||
# resulting skeleton: if you process a corner first, that corner will
|
||||
# be eroded and the skeleton will miss the arm from that corner. Pixels
|
||||
# with fewer neighbors are more "cornery" and should be processed last.
|
||||
# We use a cornerness_table lookup table where the score of a
|
||||
# configuration is the number of background (0-value) pixels in the
|
||||
# 3x3 neighbourhood
|
||||
cornerness_table = np.array([9 - np.sum(_pattern_of(index))
|
||||
for index in range(512)])
|
||||
corner_score = _table_lookup(masked_image, cornerness_table)
|
||||
|
||||
# Define arrays for inner loop
|
||||
i, j = np.mgrid[0:image.shape[0], 0:image.shape[1]]
|
||||
result = masked_image.copy()
|
||||
distance = distance[result]
|
||||
i = np.ascontiguousarray(i[result], np.int32)
|
||||
j = np.ascontiguousarray(j[result], np.int32)
|
||||
result = np.ascontiguousarray(result, np.uint8)
|
||||
|
||||
# Determine the order in which pixels are processed.
|
||||
# We use a random # for tiebreaking. Assign each pixel in the image a
|
||||
# predictable, random # so that masking doesn't affect arbitrary choices
|
||||
# of skeletons
|
||||
#
|
||||
generator = np.random.RandomState(0)
|
||||
tiebreaker = generator.permutation(np.arange(masked_image.sum()))
|
||||
order = np.lexsort((tiebreaker,
|
||||
corner_score[masked_image],
|
||||
distance))
|
||||
order = np.ascontiguousarray(order, np.int32)
|
||||
|
||||
table = np.ascontiguousarray(table, np.uint8)
|
||||
# Remove pixels not belonging to the medial axis
|
||||
_skeletonize_loop(result, i, j, order, table)
|
||||
|
||||
result = result.astype(bool)
|
||||
if not mask is None:
|
||||
result[~mask] = image[~mask]
|
||||
if return_distance:
|
||||
return result, store_distance
|
||||
else:
|
||||
return result
|
||||
|
||||
def _pattern_of(index):
|
||||
"""
|
||||
Return the pattern represented by an index value
|
||||
Byte decomposition of index
|
||||
"""
|
||||
return np.array([[index & 2**0, index & 2**1, index & 2**2],
|
||||
[index & 2**3, index & 2**4, index & 2**5],
|
||||
[index & 2**6, index & 2**7, index & 2**8]], bool)
|
||||
|
||||
|
||||
def _table_lookup(image, table):
|
||||
"""
|
||||
Perform a morphological transform on an image, directed by its
|
||||
neighbors
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
A binary image
|
||||
table : ndarray
|
||||
A 512-element table giving the transform of each pixel given
|
||||
the values of that pixel and its 8-connected neighbors.
|
||||
border_value : bool
|
||||
The value of pixels beyond the border of the image.
|
||||
|
||||
Returns
|
||||
-------
|
||||
result : ndarray of same shape as `image`
|
||||
Transformed image
|
||||
|
||||
Notes
|
||||
-----
|
||||
The pixels are numbered like this::
|
||||
|
||||
0 1 2
|
||||
3 4 5
|
||||
6 7 8
|
||||
|
||||
The index at a pixel is the sum of 2**<pixel-number> for pixels
|
||||
that evaluate to true.
|
||||
"""
|
||||
#
|
||||
# We accumulate into the indexer to get the index into the table
|
||||
# at each point in the image
|
||||
#
|
||||
if image.shape[0] < 3 or image.shape[1] < 3:
|
||||
image = image.astype(bool)
|
||||
indexer = np.zeros(image.shape, int)
|
||||
indexer[1:, 1:] += image[:-1, :-1] * 2**0
|
||||
indexer[1:, :] += image[:-1, :] * 2**1
|
||||
indexer[1:, :-1] += image[:-1, 1:] * 2**2
|
||||
|
||||
indexer[:, 1:] += image[:, :-1] * 2**3
|
||||
indexer[:, :] += image[:, :] * 2**4
|
||||
indexer[:, :-1] += image[:, 1:] * 2**5
|
||||
|
||||
indexer[:-1, 1:] += image[1:, :-1] * 2**6
|
||||
indexer[:-1, :] += image[1:, :] * 2**7
|
||||
indexer[:-1, :-1] += image[1:, 1:] * 2**8
|
||||
else:
|
||||
indexer = _table_lookup_index(np.ascontiguousarray(image, np.uint8))
|
||||
image = table[indexer]
|
||||
return image
|
||||
|
||||
|
||||
@@ -1,5 +1,5 @@
|
||||
import numpy as np
|
||||
from skimage.morphology import skeletonize
|
||||
from skimage.morphology import skeletonize, medial_axis
|
||||
import numpy.testing
|
||||
from skimage.draw import draw
|
||||
from scipy.ndimage import correlate
|
||||
@@ -90,6 +90,64 @@ class TestSkeletonize():
|
||||
blocks = correlate(result, mask, mode='constant')
|
||||
assert not numpy.any(blocks == 4)
|
||||
|
||||
class TestMedialAxis():
|
||||
def test_00_00_zeros(self):
|
||||
'''Test skeletonize on an array of all zeros'''
|
||||
result = medial_axis(np.zeros((10, 10), bool))
|
||||
assert np.all(result == False)
|
||||
|
||||
def test_00_01_zeros_masked(self):
|
||||
'''Test skeletonize on an array that is completely masked'''
|
||||
result = medial_axis(np.zeros((10, 10), bool),
|
||||
np.zeros((10, 10), bool))
|
||||
assert np.all(result == False)
|
||||
|
||||
def test_01_01_rectangle(self):
|
||||
'''Test skeletonize on a rectangle'''
|
||||
image = np.zeros((9, 15), bool)
|
||||
image[1:-1, 1:-1] = True
|
||||
#
|
||||
# The result should be four diagonals from the
|
||||
# corners, meeting in a horizontal line
|
||||
#
|
||||
expected = np.array([[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
|
||||
[0,1,0,0,0,0,0,0,0,0,0,0,0,1,0],
|
||||
[0,0,1,0,0,0,0,0,0,0,0,0,1,0,0],
|
||||
[0,0,0,1,0,0,0,0,0,0,0,1,0,0,0],
|
||||
[0,0,0,0,1,1,1,1,1,1,1,0,0,0,0],
|
||||
[0,0,0,1,0,0,0,0,0,0,0,1,0,0,0],
|
||||
[0,0,1,0,0,0,0,0,0,0,0,0,1,0,0],
|
||||
[0,1,0,0,0,0,0,0,0,0,0,0,0,1,0],
|
||||
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]], bool)
|
||||
result = medial_axis(image)
|
||||
assert np.all(result == expected)
|
||||
result, distance = medial_axis(image, return_distance=True)
|
||||
assert distance.max() == 4
|
||||
|
||||
def test_01_02_hole(self):
|
||||
'''Test skeletonize on a rectangle with a hole in the middle'''
|
||||
image = np.zeros((9, 15), bool)
|
||||
image[1:-1, 1:-1] = True
|
||||
image[4, 4:-4] = False
|
||||
expected = np.array([[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
|
||||
[0,1,0,0,0,0,0,0,0,0,0,0,0,1,0],
|
||||
[0,0,1,1,1,1,1,1,1,1,1,1,1,0,0],
|
||||
[0,0,1,0,0,0,0,0,0,0,0,0,1,0,0],
|
||||
[0,0,1,0,0,0,0,0,0,0,0,0,1,0,0],
|
||||
[0,0,1,0,0,0,0,0,0,0,0,0,1,0,0],
|
||||
[0,0,1,1,1,1,1,1,1,1,1,1,1,0,0],
|
||||
[0,1,0,0,0,0,0,0,0,0,0,0,0,1,0],
|
||||
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]],bool)
|
||||
result = medial_axis(image)
|
||||
assert np.all(result == expected)
|
||||
|
||||
def test_narrow_image(self):
|
||||
"""Test skeletonize on a 1-pixel thin strip"""
|
||||
image = np.zeros((1, 5), bool)
|
||||
image[:, 1:-1] = True
|
||||
result = medial_axis(image)
|
||||
assert np.all(result == image)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
np.testing.run_module_suite()
|
||||
|
||||
@@ -23,6 +23,7 @@ def configuration(parent_package='', top_path=None):
|
||||
config.add_extension('_project', sources=['_project.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
|
||||
return config
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
Reference in New Issue
Block a user