mirror of
https://github.com/wassname/scikit-image.git
synced 2026-08-16 11:27:48 +08:00
Juan Nunez-Iglesias
@@ -221,3 +221,6 @@
|
||||
|
||||
- Egor Panfilov
|
||||
Inpainting with biharmonic equation
|
||||
|
||||
- Evgeni Burovski
|
||||
Adaptation of ImageJ 3D skeletonization algorithm.
|
||||
|
||||
@@ -119,6 +119,9 @@ Library:
|
||||
Extension: skimage.morphology._skeletonize_cy
|
||||
Sources:
|
||||
skimage/morphology/_skeletonize_cy.pyx
|
||||
Extension: skimage.morphology._skeletonize_3d_cy
|
||||
Sources:
|
||||
skimage/morphology/_skeletonize_3d_cy.pyx
|
||||
Extension: skimage.transform._radon_transform
|
||||
Sources:
|
||||
skimage/transform/_radon_transform.pyx
|
||||
|
||||
+1
-1
@@ -91,6 +91,6 @@ def process_tempita_pyx(fromfile):
|
||||
if not fromfile.endswith('.pyx.in'):
|
||||
raise ValueError("Unexpected extension of %s." % fromfile)
|
||||
|
||||
pyxfile = os.path.splitext(fromfile)[0] + '.pyx'
|
||||
pyxfile = os.path.splitext(fromfile)[0] # split off the .in ending
|
||||
with open(pyxfile, "w") as f:
|
||||
f.write(pyxcontent)
|
||||
|
||||
Binary file not shown.
@@ -6,6 +6,7 @@ from .selem import (square, rectangle, diamond, disk, cube, octahedron, ball,
|
||||
octagon, star)
|
||||
from .watershed import watershed
|
||||
from ._skeletonize import skeletonize, medial_axis
|
||||
from ._skeletonize_3d import skeletonize_3d
|
||||
from .convex_hull import convex_hull_image, convex_hull_object
|
||||
from .greyreconstruct import reconstruction
|
||||
from .misc import remove_small_objects, remove_small_holes
|
||||
@@ -34,6 +35,7 @@ __all__ = ['binary_erosion',
|
||||
'label',
|
||||
'watershed',
|
||||
'skeletonize',
|
||||
'skeletonize_3d',
|
||||
'medial_axis',
|
||||
'convex_hull_image',
|
||||
'convex_hull_object',
|
||||
|
||||
@@ -0,0 +1,77 @@
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
import numpy as np
|
||||
from ..util import img_as_ubyte, crop
|
||||
from ._skeletonize_3d_cy import _compute_thin_image
|
||||
|
||||
|
||||
def skeletonize_3d(img):
|
||||
"""Compute the skeleton of a binary image.
|
||||
|
||||
Thinning is used to reduce each connected component in a binary image
|
||||
to a single-pixel wide skeleton.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
img : ndarray, 2D or 3D
|
||||
A binary image containing the objects to be skeletonized. Zeros
|
||||
represent background, nonzero values are foreground.
|
||||
|
||||
Returns
|
||||
-------
|
||||
skeleton : ndarray
|
||||
The thinned image.
|
||||
|
||||
See also
|
||||
--------
|
||||
skeletonize, medial_axis
|
||||
|
||||
Notes
|
||||
-----
|
||||
The method of [Lee94]_ uses an octree data structure to examine a 3x3x3
|
||||
neighborhood of a pixel. The algorithm proceeds by iteratively sweeping
|
||||
over the image, and removing pixels at each iteration until the image
|
||||
stops changing. Each iteration consists of two steps: first, a list of
|
||||
candidates for removal is assembled; then pixels from this list are
|
||||
rechecked sequentially, to better preserve connectivity of the image.
|
||||
|
||||
The algorithm this function implements is different from the algorithms
|
||||
used by either `skeletonize` or `medial_axis`, thus for 2D images the
|
||||
results produced by this function are generally different.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [Lee94] T.-C. Lee, R.L. Kashyap and C.-N. Chu, Building skeleton models
|
||||
via 3-D medial surface/axis thinning algorithms.
|
||||
Computer Vision, Graphics, and Image Processing, 56(6):462-478, 1994.
|
||||
|
||||
"""
|
||||
# make sure the image is 3D or 2D
|
||||
if img.ndim < 2 or img.ndim > 3:
|
||||
raise ValueError("skeletonize_3d can only handle 2D or 3D images; "
|
||||
"got img.ndim = %s instead." % img.ndim)
|
||||
|
||||
img = np.ascontiguousarray(img)
|
||||
img = img_as_ubyte(img, force_copy=False)
|
||||
|
||||
# make an in image 3D and pad it w/ zeros to simplify dealing w/ boundaries
|
||||
# NB: careful here to not clobber the original *and* minimize copying
|
||||
img_o = img
|
||||
if img.ndim == 2:
|
||||
img_o = img[np.newaxis, ...]
|
||||
img_o = np.pad(img_o, pad_width=1, mode='constant')
|
||||
|
||||
# normalize to binary
|
||||
maxval = img_o.max()
|
||||
img_o[img_o != 0] = 1
|
||||
|
||||
# do the computation
|
||||
img_o = np.asarray(_compute_thin_image(img_o))
|
||||
|
||||
# crop it back and restore the original intensity range
|
||||
img_o = crop(img_o, crop_width=1)
|
||||
if img.ndim == 2:
|
||||
img_o = img_o[0]
|
||||
img_o *= maxval
|
||||
|
||||
return img_o
|
||||
@@ -0,0 +1,428 @@
|
||||
"""
|
||||
This is an implementation of the 2D/3D thinning algorithm
|
||||
of [Lee94]_ of binary images, based on [IAC15]_.
|
||||
|
||||
The original Java code [IAC15]_ carries the following message:
|
||||
|
||||
* This work is an implementation by Ignacio Arganda-Carreras of the
|
||||
* 3D thinning algorithm from Lee et al. "Building skeleton models via 3-D
|
||||
* medial surface/axis thinning algorithms. Computer Vision, Graphics, and
|
||||
* Image Processing, 56(6):462-478, 1994." Based on the ITK version from
|
||||
* Hanno Homann <a href="http://hdl.handle.net/1926/1292"> http://hdl.handle.net/1926/1292</a>
|
||||
* <p>
|
||||
* More information at Skeletonize3D homepage:
|
||||
* http://fiji.sc/Skeletonize3D
|
||||
*
|
||||
* @version 1.0 11/13/2015 (unique BSD licensed version for scikit-image)
|
||||
* @author Ignacio Arganda-Carreras (iargandacarreras at gmail.com)
|
||||
|
||||
References
|
||||
----------
|
||||
.. [Lee94] T.-C. Lee, R.L. Kashyap and C.-N. Chu, Building skeleton models
|
||||
via 3-D medial surface/axis thinning algorithms.
|
||||
Computer Vision, Graphics, and Image Processing, 56(6):462-478, 1994.
|
||||
|
||||
.. [IAC15] Ignacio Arganda-Carreras, 2015. Skeletonize3D plugin for ImageJ(C).
|
||||
http://fiji.sc/Skeletonize3D
|
||||
|
||||
"""
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
from libc.string cimport memcpy
|
||||
|
||||
import numpy as np
|
||||
from numpy cimport npy_intp, npy_uint8, ndarray
|
||||
cimport cython
|
||||
|
||||
ctypedef npy_uint8 pixel_type
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
def _compute_thin_image(pixel_type[:, :, ::1] img not None):
|
||||
"""Compute a thin image.
|
||||
|
||||
Loop through the image multiple times, removing "simple" points, i.e.
|
||||
those point which can be removed without changing local connectivity in the
|
||||
3x3x3 neighborhood of a point.
|
||||
|
||||
This routine implements the two-pass algorithm of [Lee94]_. Namely,
|
||||
for each of the six border types (positive and negative x-, y- and z-),
|
||||
the algorithm first collects all possibly deletable points, and then
|
||||
performs a sequential rechecking.
|
||||
|
||||
The input, `img`, is assumed to be a 3D binary image in the
|
||||
(p, r, c) format [i.e., C ordered array], filled by zeros (background) and
|
||||
ones. Furthermore, `img` is assumed to be padded by zeros from all
|
||||
directions --- this way the zero boundary conditions are automatic
|
||||
and there is need to guard against out-of-bounds access.
|
||||
|
||||
"""
|
||||
cdef:
|
||||
int unchanged_borders = 0, curr_border, num_borders
|
||||
int borders[6]
|
||||
npy_intp p, r, c
|
||||
bint no_change
|
||||
|
||||
list simple_border_points
|
||||
Py_ssize_t num_border_points, i, j
|
||||
(npy_intp, npy_intp, npy_intp) point
|
||||
|
||||
pixel_type neighb[27]
|
||||
# loop over the six directions in this order (for consistency with ImageJ)
|
||||
borders[:] = [4, 3, 2, 1, 5, 6]
|
||||
|
||||
# no need to worry about the z direction if the original image is 2D.
|
||||
if img.shape[0] == 3:
|
||||
num_borders = 4
|
||||
else:
|
||||
num_borders = 6
|
||||
|
||||
# loop through the image several times until there is no change for all
|
||||
# the six border types
|
||||
while unchanged_borders < num_borders:
|
||||
unchanged_borders = 0
|
||||
for j in range(num_borders):
|
||||
curr_border = borders[j]
|
||||
|
||||
simple_border_points = find_simple_point_candidates(img, curr_border)
|
||||
|
||||
# sequential re-checking to preserve connectivity when deleting
|
||||
# in a parallel way
|
||||
no_change = True
|
||||
num_border_points = len(simple_border_points)
|
||||
for i in range(num_border_points):
|
||||
point = simple_border_points[i]
|
||||
p, r, c = point
|
||||
get_neighborhood(img, p, r, c, neighb)
|
||||
if is_simple_point(neighb):
|
||||
img[p, r, c] = 0
|
||||
no_change = False
|
||||
|
||||
if no_change:
|
||||
unchanged_borders += 1
|
||||
|
||||
return np.asarray(img)
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
cdef list find_simple_point_candidates(pixel_type[:, :, ::1] img,
|
||||
int curr_border):
|
||||
"""Inner loop of compute_thin_image.
|
||||
|
||||
The algorithm of [Lee94]_ proceeds in two steps: (1) six directions are
|
||||
checked for simple border points to remove, and (2) these candidates are
|
||||
sequentially rechecked, see Sec 3 of [Lee94]_ for rationale and discussion.
|
||||
|
||||
This routine implements the first step above: it loops over the image
|
||||
for a given direction and assembles candidates for removal.
|
||||
|
||||
"""
|
||||
# This routine looks like it could be nogil, but actually it cannot be,
|
||||
# because of `simple_border_points` being a python list which is being
|
||||
# mutated.
|
||||
cdef:
|
||||
list simple_border_points = []
|
||||
pixel_type neighborhood[27]
|
||||
npy_intp p, r, c
|
||||
bint is_border_pt
|
||||
(npy_intp, npy_intp, npy_intp) point
|
||||
|
||||
# rebind a global name to avoid lookup. The table is filled in
|
||||
# at import time.
|
||||
int[::1] Euler_LUT = LUT
|
||||
|
||||
# loop through the image
|
||||
# NB: each loop is from 1 to size-1: img is padded from all sides
|
||||
for p in range(1, img.shape[0] - 1):
|
||||
for r in range(1, img.shape[1] - 1):
|
||||
for c in range(1, img.shape[2] - 1):
|
||||
|
||||
# check if pixel is foreground
|
||||
if img[p, r, c] != 1:
|
||||
continue
|
||||
|
||||
is_border_pt = (curr_border == 1 and img[p, r, c-1] == 0 or #N
|
||||
curr_border == 2 and img[p, r, c+1] == 0 or #S
|
||||
curr_border == 3 and img[p, r+1, c] == 0 or #E
|
||||
curr_border == 4 and img[p, r-1, c] == 0 or #W
|
||||
curr_border == 5 and img[p+1, r, c] == 0 or #U
|
||||
curr_border == 6 and img[p-1, r, c] == 0) #B
|
||||
if not is_border_pt:
|
||||
# current point is not deletable
|
||||
continue
|
||||
|
||||
get_neighborhood(img, p, r, c, neighborhood)
|
||||
|
||||
# check if (p, r, c) can be deleted:
|
||||
# * it must not be an endpoint;
|
||||
# * it must be Euler invariant (condition 1 in [Lee94]_); and
|
||||
# * it must be simple (i.e., its deletion does not change
|
||||
# connectivity in the 3x3x3 neighborhood)
|
||||
# this is conditions 2 and 3 in [Lee94]_
|
||||
if (is_endpoint(neighborhood) or
|
||||
not is_Euler_invariant(neighborhood, Euler_LUT) or
|
||||
not is_simple_point(neighborhood)):
|
||||
continue
|
||||
|
||||
# ok, add (p, r, c) to the list of simple border points
|
||||
point = p, r, c
|
||||
simple_border_points.append(point)
|
||||
return simple_border_points
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
cdef void get_neighborhood(pixel_type[:, :, ::1] img,
|
||||
npy_intp p, npy_intp r, npy_intp c,
|
||||
pixel_type neighborhood[]):
|
||||
"""Get the neighborhood of a pixel.
|
||||
|
||||
Assume zero boundary conditions.
|
||||
Image is already padded, so no out-of-bounds checking.
|
||||
|
||||
For the numbering of points see Fig. 1a. of [Lee94]_, where the numbers
|
||||
do *not* include the center point itself. OTOH, this numbering below
|
||||
includes it as number 13. The latter is consistent with [IAC15]_.
|
||||
"""
|
||||
neighborhood[0] = img[p-1, r-1, c-1]
|
||||
neighborhood[1] = img[p-1, r, c-1]
|
||||
neighborhood[2] = img[p-1, r+1, c-1]
|
||||
|
||||
neighborhood[ 3] = img[p-1, r-1, c]
|
||||
neighborhood[ 4] = img[p-1, r, c]
|
||||
neighborhood[ 5] = img[p-1, r+1, c]
|
||||
|
||||
neighborhood[ 6] = img[p-1, r-1, c+1]
|
||||
neighborhood[ 7] = img[p-1, r, c+1]
|
||||
neighborhood[ 8] = img[p-1, r+1, c+1]
|
||||
|
||||
neighborhood[ 9] = img[p, r-1, c-1]
|
||||
neighborhood[10] = img[p, r, c-1]
|
||||
neighborhood[11] = img[p, r+1, c-1]
|
||||
|
||||
neighborhood[12] = img[p, r-1, c]
|
||||
neighborhood[13] = img[p, r, c]
|
||||
neighborhood[14] = img[p, r+1, c]
|
||||
|
||||
neighborhood[15] = img[p, r-1, c+1]
|
||||
neighborhood[16] = img[p, r, c+1]
|
||||
neighborhood[17] = img[p, r+1, c+1]
|
||||
|
||||
neighborhood[18] = img[p+1, r-1, c-1]
|
||||
neighborhood[19] = img[p+1, r, c-1]
|
||||
neighborhood[20] = img[p+1, r+1, c-1]
|
||||
|
||||
neighborhood[21] = img[p+1, r-1, c]
|
||||
neighborhood[22] = img[p+1, r, c]
|
||||
neighborhood[23] = img[p+1, r+1, c]
|
||||
|
||||
neighborhood[24] = img[p+1, r-1, c+1]
|
||||
neighborhood[25] = img[p+1, r, c+1]
|
||||
neighborhood[26] = img[p+1, r+1, c+1]
|
||||
|
||||
|
||||
###### look-up tables
|
||||
def fill_Euler_LUT():
|
||||
""" Look-up table for preserving Euler characteristic.
|
||||
|
||||
This is column $\delta G_{26}$ of Table 2 of [Lee94]_.
|
||||
"""
|
||||
cdef int arr[128]
|
||||
arr[:] = [1, -1, -1, 1, -3, -1, -1, 1, -1, 1, 1, -1, 3, 1, 1, -1, -3, -1,
|
||||
3, 1, 1, -1, 3, 1, -1, 1, 1, -1, 3, 1, 1, -1, -3, 3, -1, 1, 1,
|
||||
3, -1, 1, -1, 1, 1, -1, 3, 1, 1, -1, 1, 3, 3, 1, 5, 3, 3, 1,
|
||||
-1, 1, 1, -1, 3, 1, 1, -1, -7, -1, -1, 1, -3, -1, -1, 1, -1,
|
||||
1, 1, -1, 3, 1, 1, -1, -3, -1, 3, 1, 1, -1, 3, 1, -1, 1, 1,
|
||||
-1, 3, 1, 1, -1, -3, 3, -1, 1, 1, 3, -1, 1, -1, 1, 1, -1, 3,
|
||||
1, 1, -1, 1, 3, 3, 1, 5, 3, 3, 1, -1, 1, 1, -1, 3, 1, 1, -1]
|
||||
cdef ndarray LUT = np.zeros(256, dtype=np.intc)
|
||||
LUT[1::2] = arr
|
||||
return LUT
|
||||
cdef int[::1] LUT = fill_Euler_LUT()
|
||||
|
||||
|
||||
# Fill the look-up table for indexing octants for computing the Euler
|
||||
# characteristic. See is_Euler_invariant routine below.
|
||||
{{py:
|
||||
_neighb_idx = [[2, 1, 11, 10, 5, 4, 14], # NEB
|
||||
[0, 9, 3, 12, 1, 10, 4], # NWB
|
||||
[8, 7, 17, 16, 5, 4, 14], # SEB
|
||||
[6, 15, 7, 16, 3, 12, 4], # SWB
|
||||
[20, 23, 19, 22, 11, 14, 10], # NEU
|
||||
[18, 21, 9, 12, 19, 22, 10], # NWU
|
||||
[26, 23, 17, 14, 25, 22, 16], # SEU
|
||||
[24, 25, 15, 16, 21, 22, 12], # SWU
|
||||
]
|
||||
}}
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
cdef bint is_Euler_invariant(pixel_type neighbors[],
|
||||
int[::1] lut):
|
||||
"""Check if a point is Euler invariant.
|
||||
|
||||
Calculate Euler characteristic for each octant and sum up.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
neighbors
|
||||
neighbors of a point
|
||||
lut
|
||||
The look-up table for preserving the Euler characteristic.
|
||||
|
||||
Returns
|
||||
-------
|
||||
bool (C bool, that is)
|
||||
|
||||
"""
|
||||
cdef int n, euler_char = 0
|
||||
{{for _octant in range(8)}}
|
||||
|
||||
# octant {{_octant}}:
|
||||
n = 1
|
||||
{{for _j in range(7):}}
|
||||
{{py: _idx = _neighb_idx[_octant][_j]}}
|
||||
if neighbors[{{_idx}}] == 1:
|
||||
n |= {{1 << (7 - _j)}}
|
||||
|
||||
{{endfor}}
|
||||
euler_char += lut[n]
|
||||
{{endfor}}
|
||||
return euler_char == 0
|
||||
|
||||
|
||||
cdef inline bint is_endpoint(pixel_type neighbors[]):
|
||||
"""An endpoint has exactly one neighbor in the 26-neighborhood.
|
||||
"""
|
||||
# The center pixel is counted, thus r.h.s. is 2
|
||||
cdef int s = 0, j
|
||||
for j in range(27):
|
||||
s += neighbors[j]
|
||||
return s == 2
|
||||
|
||||
|
||||
cdef bint is_simple_point(pixel_type neighbors[]):
|
||||
"""Check is a point is a Simple Point.
|
||||
|
||||
A point is simple iff its deletion does not change connectivity in
|
||||
the 3x3x3 neighborhood. (cf conditions 2 and 3 in [Lee94]_).
|
||||
|
||||
This method is named "N(v)_labeling" in [Lee94]_.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
neighbors : uint8 C array, shape(27,)
|
||||
neighbors of the point
|
||||
|
||||
Returns
|
||||
-------
|
||||
bool
|
||||
Whether the point is simple or not.
|
||||
|
||||
"""
|
||||
# copy neighbors for labeling
|
||||
# ignore center pixel (i=13) when counting (see [Lee94]_)
|
||||
cdef pixel_type cube[26]
|
||||
memcpy(cube, neighbors, 13*sizeof(pixel_type))
|
||||
memcpy(cube+13, neighbors+14, 13*sizeof(pixel_type))
|
||||
|
||||
# set initial label
|
||||
cdef int label = 2, i
|
||||
|
||||
# for all point in the neighborhood
|
||||
for i in range(26):
|
||||
if cube[i] == 1:
|
||||
# voxel has not been labeled yet
|
||||
# start recursion with any octant that contains the point i
|
||||
if i in (0, 1, 3, 4, 9, 10, 12):
|
||||
octree_labeling(1, label, cube)
|
||||
elif i in (2, 5, 11, 13):
|
||||
octree_labeling(2, label, cube)
|
||||
elif i in (6, 7, 14, 15):
|
||||
octree_labeling(3, label, cube)
|
||||
elif i in (8, 16):
|
||||
octree_labeling(4, label, cube)
|
||||
elif i in (17, 18, 20, 21):
|
||||
octree_labeling(5, label, cube)
|
||||
elif i in (19, 22):
|
||||
octree_labeling(6, label, cube)
|
||||
elif i in (23, 24):
|
||||
octree_labeling(7, label, cube)
|
||||
elif i == 25:
|
||||
octree_labeling(8, label, cube)
|
||||
label += 1
|
||||
if label - 2 >= 2:
|
||||
return False
|
||||
return True
|
||||
|
||||
|
||||
# Octree structure for labeling in `octree_labeling` routine below.
|
||||
# NB: this is only available at build time, and is used by Tempita templating.
|
||||
{{py:
|
||||
_octree = [
|
||||
# octant 1
|
||||
([0, 1, 3, 4, 9, 10, 12],
|
||||
[[], [2], [3], [2, 3, 4], [5], [2, 5, 6], [3, 5, 7]]),
|
||||
# octant 2
|
||||
([1, 4, 10, 2, 5, 11, 13],
|
||||
[[1], [1, 3, 4], [1, 5, 6], [], [4], [6], [4, 6, 8]]),
|
||||
# octant 3
|
||||
([3, 4, 12, 6, 7, 14, 15],
|
||||
[[1], [1, 2, 4], [1, 5, 7], [], [4], [7], [4, 7, 8]]),
|
||||
# octant 4
|
||||
([4, 5, 13, 7, 15, 8, 16],
|
||||
[[1, 2, 3], [2], [2, 6, 8], [3], [3, 7, 8], [], [8]]),
|
||||
# octant 5
|
||||
([9, 10, 12, 17, 18, 20, 21],
|
||||
[[1], [1, 2, 6], [1, 3, 7], [], [6], [7], [6, 7, 8]]),
|
||||
# octant 6
|
||||
([10, 11, 13, 18, 21, 19, 22],
|
||||
[[1, 2, 5], [2], [2, 4, 8], [5], [5, 7, 8], [], [8]]),
|
||||
# octant 7
|
||||
([12, 14, 15, 20, 21, 23, 24],
|
||||
[[1, 3, 5], [3], [3, 4, 8], [5], [5, 6, 8], [], [8]]),
|
||||
# octant 8
|
||||
([13, 15, 16, 21, 22, 24, 25],
|
||||
[[2, 4, 6], [3, 4, 7], [4], [5, 6, 7], [6], [7], []])
|
||||
]
|
||||
}}
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
cdef void octree_labeling(int octant, int label, pixel_type cube[]):
|
||||
"""This is a recursive method that calculates the number of connected
|
||||
components in the 3D neighborhood after the center pixel would
|
||||
have been removed.
|
||||
|
||||
See Figs. 6 and 7 of [Lee94]_ for the values of indices.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
octant : int
|
||||
octant index
|
||||
label : int
|
||||
the current label of the center point
|
||||
cube : uint8 C array, shape(26,)
|
||||
local neighborhood of the point
|
||||
|
||||
"""
|
||||
# This routine checks if there are points in the octant with value 1
|
||||
# Then sets points in this octant to current label
|
||||
# and recursive labeling of adjacent octants.
|
||||
#
|
||||
# Below, leading underscore means build-time variables.
|
||||
{{for _oct in range(1, 9)}}
|
||||
|
||||
if octant == {{_oct}}:
|
||||
{{py: _indices, _list_octants = _octree[_oct-1]}}
|
||||
{{for _idx, _new_octants in zip(_indices, _list_octants)}}
|
||||
if cube[{{_idx}}] == 1:
|
||||
cube[{{_idx}}] = label
|
||||
{{for _new_octant in _new_octants}}
|
||||
octree_labeling({{_new_octant}}, label, cube)
|
||||
{{endfor}}
|
||||
{{endfor}}
|
||||
{{endfor}}
|
||||
@@ -16,6 +16,7 @@ def configuration(parent_package='', top_path=None):
|
||||
cython(['_skeletonize_cy.pyx'], working_path=base_path)
|
||||
cython(['_convex_hull.pyx'], working_path=base_path)
|
||||
cython(['_greyreconstruct.pyx'], working_path=base_path)
|
||||
cython(['_skeletonize_3d_cy.pyx.in'], working_path=base_path)
|
||||
|
||||
config.add_extension('_watershed', sources=['_watershed.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
@@ -25,6 +26,8 @@ def configuration(parent_package='', top_path=None):
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_greyreconstruct', sources=['_greyreconstruct.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_skeletonize_3d_cy', sources=['_skeletonize_3d_cy.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
return config
|
||||
|
||||
|
||||
@@ -0,0 +1,186 @@
|
||||
from __future__ import division, print_function, absolute_import
|
||||
|
||||
import os
|
||||
import warnings
|
||||
|
||||
import numpy as np
|
||||
from numpy.testing import (assert_equal, run_module_suite, assert_raises,
|
||||
assert_)
|
||||
|
||||
import scipy.ndimage as ndi
|
||||
|
||||
import skimage
|
||||
from skimage import io, draw, data_dir
|
||||
from skimage.data import binary_blobs
|
||||
from skimage.util import img_as_ubyte
|
||||
|
||||
from skimage.morphology import skeletonize_3d
|
||||
|
||||
|
||||
# basic behavior tests (mostly copied over from 2D skeletonize)
|
||||
|
||||
def test_skeletonize_wrong_dim():
|
||||
im = np.zeros(5, dtype=np.uint8)
|
||||
assert_raises(ValueError, skeletonize_3d, im)
|
||||
|
||||
im = np.zeros((5, 5, 5, 5), dtype=np.uint8)
|
||||
assert_raises(ValueError, skeletonize_3d, im)
|
||||
|
||||
|
||||
def test_skeletonize_1D():
|
||||
# a corner case of an image of a shape(1, N)
|
||||
im = np.ones((5, 1), dtype=np.uint8)
|
||||
res = skeletonize_3d(im)
|
||||
assert_equal(res, im)
|
||||
|
||||
|
||||
def test_skeletonize_no_foreground():
|
||||
im = np.zeros((5, 5), dtype=np.uint8)
|
||||
result = skeletonize_3d(im)
|
||||
assert_equal(result, im)
|
||||
|
||||
|
||||
def test_skeletonize_all_foreground():
|
||||
im = np.ones((3, 4), dtype=np.uint8)
|
||||
assert_equal(skeletonize_3d(im),
|
||||
np.array([[0, 0, 0, 0],
|
||||
[1, 1, 1, 1],
|
||||
[0, 0, 0, 0]], dtype=np.uint8))
|
||||
|
||||
|
||||
def test_skeletonize_single_point():
|
||||
im = np.zeros((5, 5), dtype=np.uint8)
|
||||
im[3, 3] = 1
|
||||
result = skeletonize_3d(im)
|
||||
assert_equal(result, im)
|
||||
|
||||
|
||||
def test_skeletonize_already_thinned():
|
||||
im = np.zeros((5, 5), dtype=np.uint8)
|
||||
im[3, 1:-1] = 1
|
||||
im[2, -1] = 1
|
||||
im[4, 0] = 1
|
||||
result = skeletonize_3d(im)
|
||||
assert_equal(result, im)
|
||||
|
||||
|
||||
def test_dtype_conv():
|
||||
# check that the operation does the right thing with floats etc
|
||||
# also check non-contiguous input
|
||||
img = np.random.random((16, 16))[::2, ::2]
|
||||
img[img < 0.5] = 0
|
||||
|
||||
orig = img.copy()
|
||||
|
||||
with warnings.catch_warnings():
|
||||
# UserWarning for possible precision loss, expected
|
||||
warnings.simplefilter('ignore', UserWarning)
|
||||
res = skeletonize_3d(img)
|
||||
|
||||
assert_equal(res.dtype, np.uint8)
|
||||
assert_equal(img, orig) # operation does not clobber the original
|
||||
assert_equal(res.max(),
|
||||
img_as_ubyte(img).max()) # the intensity range is preserved
|
||||
|
||||
|
||||
def test_input():
|
||||
# check that the input is not clobbered
|
||||
# for 2D and 3D images of varying dtypes
|
||||
imgs = [np.ones((8, 8), dtype=float), np.ones((4, 8, 8), dtype=float),
|
||||
np.ones((8, 8), dtype=np.uint8), np.ones((4, 8, 8), dtype=np.uint8),
|
||||
np.ones((8, 8), dtype=bool), np.ones((4, 8, 8), dtype=bool)]
|
||||
for img in imgs:
|
||||
yield check_input, img
|
||||
|
||||
|
||||
def check_input(img):
|
||||
orig = img.copy()
|
||||
with warnings.catch_warnings():
|
||||
# UserWarning for possible precision loss, expected
|
||||
warnings.simplefilter('ignore', UserWarning)
|
||||
res = skeletonize_3d(img)
|
||||
assert_equal(img, orig)
|
||||
|
||||
|
||||
def test_skeletonize_num_neighbours():
|
||||
# an empty image
|
||||
image = np.zeros((300, 300))
|
||||
|
||||
# foreground object 1
|
||||
image[10:-10, 10:100] = 1
|
||||
image[-100:-10, 10:-10] = 1
|
||||
image[10:-10, -100:-10] = 1
|
||||
|
||||
# foreground object 2
|
||||
rs, cs = draw.line(250, 150, 10, 280)
|
||||
for i in range(10):
|
||||
image[rs + i, cs] = 1
|
||||
rs, cs = draw.line(10, 150, 250, 280)
|
||||
for i in range(20):
|
||||
image[rs + i, cs] = 1
|
||||
|
||||
# foreground object 3
|
||||
ir, ic = np.indices(image.shape)
|
||||
circle1 = (ic - 135)**2 + (ir - 150)**2 < 30**2
|
||||
circle2 = (ic - 135)**2 + (ir - 150)**2 < 20**2
|
||||
image[circle1] = 1
|
||||
image[circle2] = 0
|
||||
result = skeletonize_3d(image)
|
||||
|
||||
# there should never be a 2x2 block of foreground pixels in a skeleton
|
||||
mask = np.array([[1, 1],
|
||||
[1, 1]], np.uint8)
|
||||
blocks = ndi.correlate(result, mask, mode='constant')
|
||||
assert_(not np.any(blocks == 4))
|
||||
|
||||
|
||||
def test_two_hole_image():
|
||||
# test a simple 2D image against FIJI
|
||||
img_o = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0],
|
||||
[0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
|
||||
[0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0],
|
||||
[0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0],
|
||||
[0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0],
|
||||
[0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0],
|
||||
[0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0],
|
||||
[0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
|
||||
[0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
|
||||
[0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
|
||||
[0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],
|
||||
dtype=np.uint8)
|
||||
img_f = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0],
|
||||
[0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0],
|
||||
[0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
|
||||
[0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
|
||||
[0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
|
||||
[0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
|
||||
[0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0],
|
||||
[0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],
|
||||
dtype=np.uint8)
|
||||
res = skeletonize_3d(img_o)
|
||||
assert_equal(res, img_f)
|
||||
|
||||
|
||||
def test_3d_vs_fiji():
|
||||
# generate an image with blobs and compate its skeleton to
|
||||
# the skeleton generated by FIJI
|
||||
img = binary_blobs(32, 0.05, n_dim=3, seed=1234)
|
||||
img = img[:-2, ...]
|
||||
img = img.astype(np.uint8)*255
|
||||
|
||||
img_s = skeletonize_3d(img)
|
||||
img_f = io.imread(os.path.join(data_dir, "_blobs_3d_fiji_skeleton.tif"))
|
||||
assert_equal(img_s, img_f)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
run_module_suite()
|
||||
@@ -30,7 +30,7 @@ retry () {
|
||||
}
|
||||
|
||||
# add build dependencies
|
||||
echo "cython>=0.21" >> requirements.txt
|
||||
echo "cython>=0.23.4" >> requirements.txt
|
||||
|
||||
# require networkx 1.9.1 on 2.6, as 2.6 support was dropped in 1.10
|
||||
# require matplotlib 1.4.3 on 2.6, as 2.6 support was dropped in 1.5
|
||||
|
||||
Reference in New Issue
Block a user