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

Medial axis transform (skeletonization) from cell profiler

Browse Source 
* medial_axis function in morphology.skeletonize
 * plot_medial_axis.py example for the gallery

Compared to the skeletonize algorithm in morphology.skeletonize,
medial_axis is faster because it processes all pixel in one pass.
However, the resulting skeleton has more branches, that may be unwanted
depending on the application.
This commit is contained in:
emmanuelle
2011-10-23 11:06:56 +02:00
parent 0f31034a46
commit 736802bad4
4 changed files with 183 additions and 225 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/examples/plot_medial_transform.py
+70
View File
@@ -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 less 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 scikits.image.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()
scikits/image/morphology/_cpmorphology2.pyx
+34 -119
View File
@@ -38,30 +38,49 @@ def skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2,
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
"""
Inner loop of skeletonize function
result - on input, the image to be skeletonized, on output the skeletonized
image
i,j - the coordinates of each foreground pixel in the image
order - the index of each pixel, in the order of processing
table - the 512-element lookup table of values after transformation
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
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
@@ -74,13 +93,14 @@ def skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2,
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]:
if jj > 0 and result[ii + 1, jj - 1]:
accumulator += 64
if result[ii+1,jj]:
if result[ii + 1, jj]:
accumulator += 128
if jj < result.shape[1]-1 and result[ii+1,jj+1]:
if jj < result.shape[1] - 1 and result[ii + 1, jj + 1]:
accumulator += 256
result[ii,jj] = table[accumulator]
# 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,
@@ -198,112 +218,7 @@ def table_lookup_index(np.ndarray[dtype=np.uint8_t, ndim=2,
indexer[i + 1, j_shape - 1] += 2
return indexer
@cython.boundscheck(False)
def index_lookup(np.ndarray[dtype=np.int32_t, ndim=1,
negative_indices=False] index_i,
np.ndarray[dtype=np.int32_t, ndim=1,
negative_indices=False] index_j,
np.ndarray[dtype=np.uint32_t, ndim=2,
negative_indices=False] image,
table_in,
iterations=None):
"""
Perform a table lookup for only the indexed pixels
For morphological operations that only convert 1 to 0, the set of
resulting pixels is always a subset of the input set. Therefore, when
repeating, it will be faster to operate only on the subsets especially
when the results are 1-d or 0-d objects.
This function returns a new index_i and index_j array of the pixels
that survive the operation. The image is modified in-place to remove
the pixels that did not survive.
index_i - an array of row indexes into the image.
index_j - a similarly-shaped array of column indexes.
image - the binary image: *NOTE* add a row and column of border values
to the original image to account for pixels on the edge of the
image.
iterations - # of iterations to do, default is "forever"
The idea of index_lookup was taken from
http://blogs.mathworks.com/steve/2008/06/13/performance-optimization-for-applylut/
which, apparently, is how Matlab achieved its bwmorph speedup.
"""
cdef:
np.ndarray[dtype=np.uint8_t, ndim=1,
negative_indices=False] table = table_in.astype(np.uint8)
np.uint32_t center, hit_count, idx, indexer
np.int32_t idxi, idxj
if iterations == None:
# Worst case - remove one per iteration
iterations = len(index_i)
for i in range(iterations):
hit_count = len(index_i)
with nogil:
#
# For integer images (i.e., labels), a neighbor point is
# "background" if it doesn't match the central value. This
# lets adjacent labeled objects shrink independently of each
# other.
#
for 0 <= idx < hit_count:
idxi, idxj = index_i[idx], index_j[idx]
center = image[idxi, idxj]
indexer = ((image[idxi - 1, idxj - 1] == center) * 1 +
(image[idxi - 1, idxj] == center) * 2 +
(image[idxi - 1, idxj + 1] == center) * 4 +
(image[idxi, idxj - 1] == center) * 8 +
16 +
(image[idxi, idxj + 1] == center) * 32 +
(image[idxi + 1, idxj - 1] == center) * 64 +
(image[idxi + 1, idxj] == center) * 128 +
(image[idxi + 1, idxj + 1] == center) * 256)
if table[indexer] == 0:
# mark for deletion
index_i[idx] = -index_i[idx]
# remove marked pixels
for 0 <= idx < hit_count:
idxi, idxj = index_i[idx], index_j[idx]
if idxi < 0:
image[-idxi, idxj] = 0
index_j = index_j[index_i >= 0]
index_i = index_i[index_i >= 0]
if len(index_i) == hit_count:
break
return (index_i, index_j)
def prepare_for_index_lookup(image, border_value):
"""
Return the index arrays of "1" pixels and an image with an added border
The routine, index_lookup takes an array of i indexes, an array of
j indexes and an image guaranteed to be indexed successfully by
index_<i,j>[:] +/- 1. This routine constructs an image with added border
pixels... evilly, the index, 0 - 1, lands on the border because of Python's
negative indexing convention.
"""
if np.issubdtype(image.dtype, float):
image = image.astype(bool)
image_i, image_j = np.argwhere(image.astype(bool)).transpose().\
astype(np.int32) + 1
output_image = (np.ones(np.array(image.shape) + 2, image.dtype) \
if border_value
else np.zeros(np.array(image.shape) + 2, image.dtype))
output_image[1:image.shape[0] + 1, 1:image.shape[1] + 1] = image
return (image_i, image_j, output_image.astype(np.uint32))
def extract_from_image_lookup(orig_image, index_i, index_j):
"""
Extract only one pixel
"""
output = np.zeros(orig_image.shape, orig_image.dtype)
output[index_i - 1, index_j - 1] = orig_image[index_i - 1, index_j - 1]
return output
skimage/morphology/skeletonize.py
+70 -106
View File
@@ -1,14 +1,11 @@
"""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 import ndimage
from _cpmorphology2 import skeletonize_loop, table_lookup_index
from _cpmorphology2 import extract_from_image_lookup, \
prepare_for_index_lookup, index_lookup
# --------- Skeletonization by morphological thinning ---------
@@ -168,7 +165,7 @@ def medial_axis(image, mask=None, return_distance=False):
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
return_distance: bool, optional
If true, the distance transform is returned as well as the skeleton.
Returns
@@ -188,11 +185,24 @@ def medial_axis(image, mask=None, return_distance=False):
Notes
-----
This algorithm computes the medial axis transform of an image
as the ridges of its distance transform. First, the distance transform
is computed, then the foreground (value of 1) points are ordered by
the distance transform. In order to reduce the image to its skeleton,
a point is removed if it has more than one neighbor and if removing it
does not change the Euler number (the connectivity).
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 not a pixel according to the lookup table. Because of the
ordering, it is possible to process all pixels in only one pass.
Examples
--------
@@ -223,48 +233,60 @@ def medial_axis(image, mask=None, return_distance=False):
masked_image = image.astype(bool).copy()
masked_image[~mask] = False
#
# Lookup table - start with only positive pixels.
# Keep if # pixels in neighborhood is 2 or less
# Keep if removing the pixel results in a different connectivity
# table is independent of image
table = (_make_table(True,
np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]], bool),
np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]], bool)) &
# 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 &
(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)]) |
np.array([np.sum(_pattern_of(index)) < 3 for index in range(512)])))
# Build distance transform
distance = ndimage.distance_transform_edt(masked_image)
if return_distance:
store_distance = distance.copy()
#
# The processing order along the edge is critical to the shape of the
# 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, False, 1)
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
#
# Why fix the seed? Should we pass a random number generator instead?
np.random.seed(0)
tiebreaker = np.random.permutation(np.arange(masked_image.sum()))
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)
@@ -282,12 +304,12 @@ 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)
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, border_value, iterations = None):
def _table_lookup(image, table):
"""
Perform a morphological transform on an image, directed by its
neighbors
@@ -316,85 +338,27 @@ def _table_lookup(image, table, border_value, iterations = None):
that evaluate to true.
"""
#
# Test for a table that never transforms a zero into a one:
# We accumulate into the indexer to get the index into the table
# at each point in the image
#
center_is_zero = np.array([(x & 2**4) == 0 for x in range(2**9)])
use_index_trick = False
if (not np.any(table[center_is_zero]) and
(np.issubdtype(image.dtype, bool) or np.issubdtype(image.dtype, int))):
# Use the index trick
use_index_trick = True
invert = False
elif (np.all(table[~center_is_zero]) and np.issubdtype(image.dtype, bool)):
# All ones stay ones, invert the table and the image and do the trick
use_index_trick = True
invert = True
image = ~ image
# table index 0 -> 511 and the output is reversed
table = ~ table[511-np.arange(512)]
border_value = not border_value
if use_index_trick:
orig_image = image
index_i, index_j, image = prepare_for_index_lookup(image, border_value)
index_i, index_j = index_lookup(index_i, index_j,
image, table, iterations)
image = extract_from_image_lookup(orig_image, index_i, index_j)
if invert:
image = ~ image
return image
print(use_index_trick)
counter = 0
while counter != iterations:
counter += 1
#
# 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
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, 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))
if border_value:
indexer[0,:] |= 2**0 + 2**1 + 2**2
indexer[-1,:] |= 2**6 + 2**7 + 2**8
indexer[:,0] |= 2**0 + 2**3 + 2**6
indexer[:,-1] |= 2**2 + 2**5 + 2**8
new_image = table[indexer]
if np.all(new_image == image):
break
image = new_image
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
def _make_table(value, pattern, care=np.ones((3,3),bool)):
'''Return a table suitable for table_lookup
value - set all table entries matching "pattern" to "value", all others
to not "value"
pattern - a 3x3 boolean array with the pattern to match
care - a 3x3 boolean array where each value is true if the pattern
must match at that position and false if we don't care if
the pattern matches at that position.
'''
def fn(index, p, i, j):
'''Return true if bit position "p" in index matches pattern'''
return ((((index & 2**p) > 0) == pattern[i, j]) or not care[i, j])
return np.array([value
if (fn(i, 0, 0, 0) and fn(i, 1, 0, 1) and fn(i, 2, 0, 2)
and fn(i, 3, 1, 0) and fn(i, 4, 1, 1) and fn(i, 5, 1, 2)
and fn(i, 6, 2, 0) and fn(i, 7, 2, 1) and fn(i, 8, 2, 2))
else not value
for i in range(512)], bool)
skimage/morphology/tests/test_skeletonize.py
+9
View File
@@ -121,6 +121,8 @@ class TestMedialAxis():
[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'''
@@ -139,6 +141,13 @@ class TestMedialAxis():
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()