mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-08 22:23:07 +08:00
First integration of cell profiler medial axis transform
(skeletonization in cell profiler)
This commit is contained in:
emmanuelle
@@ -0,0 +1,309 @@
|
||||
'''_cpmorphology2.pyx - support routines for cpmorphology in Cython
|
||||
|
||||
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
|
||||
|
||||
cdef extern from "Python.h":
|
||||
ctypedef int Py_intptr_t
|
||||
|
||||
cdef extern from "numpy/arrayobject.h":
|
||||
ctypedef class numpy.ndarray [object PyArrayObject]:
|
||||
cdef char *data
|
||||
cdef Py_intptr_t *dimensions
|
||||
cdef Py_intptr_t *strides
|
||||
cdef void import_array()
|
||||
cdef int PyArray_ITEMSIZE(np.ndarray)
|
||||
|
||||
import_array()
|
||||
|
||||
@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
|
||||
|
||||
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
|
||||
|
||||
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.
|
||||
'''
|
||||
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]
|
||||
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
|
||||
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
|
||||
|
||||
@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
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -4,7 +4,13 @@ objects in an image.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
from scipy.ndimage import correlate
|
||||
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 ---------
|
||||
|
||||
def skeletonize(image):
|
||||
"""Return the skeleton of a binary image.
|
||||
@@ -24,6 +30,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 +117,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 +136,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 +149,252 @@ 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. 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).
|
||||
|
||||
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
|
||||
#
|
||||
# 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)) &
|
||||
(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)])))
|
||||
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
|
||||
# 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.
|
||||
#
|
||||
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)
|
||||
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)
|
||||
#
|
||||
# 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()))
|
||||
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, border_value, iterations = None):
|
||||
"""
|
||||
Perform a morphological transform on an image, directed by its
|
||||
neighbors
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image - a binary image
|
||||
table - a 512-element table giving the transform of each pixel given
|
||||
the values of that pixel and its 8-connected neighbors.
|
||||
border_value - the value of pixels beyond the border of the image.
|
||||
This should test as True or False.
|
||||
|
||||
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.
|
||||
"""
|
||||
#
|
||||
# Test for a table that never transforms a zero into a one:
|
||||
#
|
||||
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
|
||||
|
||||
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
|
||||
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)
|
||||
|
||||
|
||||
|
||||
@@ -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,55 @@ 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)
|
||||
|
||||
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)
|
||||
|
||||
|
||||
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