From 88a9bcc874bfe8ce7dc1a93653236bf17592a3c7 Mon Sep 17 00:00:00 2001 From: Evgeni Burovski Date: Tue, 2 Feb 2016 15:52:39 +0000 Subject: [PATCH] DOC: skel3d: add basic docstrings --- skimage/morphology/_skel.pyx | 81 ++++++++++++++++++++++++++++++++---- skimage/morphology/skel.py | 3 ++ 2 files changed, 77 insertions(+), 7 deletions(-) diff --git a/skimage/morphology/_skel.pyx b/skimage/morphology/_skel.pyx index b598f3db..937b8a8e 100644 --- a/skimage/morphology/_skel.pyx +++ b/skimage/morphology/_skel.pyx @@ -1,3 +1,34 @@ +""" +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 http://hdl.handle.net/1926/1292 + *

+ * 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) + +Porting to Cython was done by Evgeni Burovski (evgeny.burovskiy@gmail.com). + +References +---------- + +.. [Lee94] 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 + +.. [IAC15] Ignacio Arganda-Carreras, 2015. Skeletonize3D plugin for ImageJ(C). + http://fiji.sc/Skeletonize3D + +""" from __future__ import division, print_function, absolute_import import numpy as np @@ -14,8 +45,12 @@ cdef void get_neighborhood(pixel_type[:, :, ::1] img, pixel_type[::1] neighborhood): """Get the neighborhood of a pixel. - Assume zero boundary conditions. Image is already padded, so no - out-of-bounds checking. + 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] @@ -63,6 +98,10 @@ NUMPOINTS_LUT = fill_numpoints_LUT() def fill_Euler_LUT(): + """ Look-up table for preserving Euler characteristic. + + This is column $\delta G_{26}$ of Table 2 of [Lee94]. + """ LUT = np.zeros(256, dtype=np.intc) LUT[1] = 1 @@ -217,6 +256,7 @@ _neib_idx[SEU, ...] = [26, 23, 17, 14, 25, 22, 16] _neib_idx[SWU, ...] = [24, 25, 15, 16, 21, 22, 12] cdef int[:, ::1] neib_idx = _neib_idx + @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) @@ -266,7 +306,7 @@ cdef bint is_Euler_invariant(pixel_type[::1] neighbors): Returns ------- - bool + bool (C bool, that is) """ cdef int octant, n @@ -282,7 +322,7 @@ cdef bint is_Euler_invariant(pixel_type[::1] neighbors): cdef bint is_simple_point(pixel_type[::1] neighbors): """Check is a point is a Simple Point. - This method is named 'N(v)_labeling' in [Lee94]. + This method is named "N(v)_labeling" in [Lee94]. Outputs the number of connected objects in a neighborhood of a point after this point would have been removed. @@ -344,6 +384,8 @@ cdef void octree_labeling(int octant, int label, pixel_type[::1] cube): 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 @@ -588,8 +630,17 @@ cdef list _loop_through(pixel_type[:, :, ::1] img, int curr_border): """Inner loop of compute_thin_image. - return simple_border_points as a list to be rechecked sequentially. + 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[::1] neighborhood = np.zeros(27, dtype=np.uint8) @@ -618,7 +669,7 @@ cdef list _loop_through(pixel_type[:, :, ::1] img, get_neighborhood(img, p, r, c, neighborhood) - # check if (p, r, c) is an endpoint: endpoints are not deletable. + # check if (p, r, c) is an endpoint (then it's not deletable.) if is_endpoint(neighborhood): continue @@ -641,7 +692,24 @@ cdef list _loop_through(pixel_type[:, :, ::1] img, @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 algorthim 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 authomatic + and there is need to guard against out-of-bounds access. + + """ cdef: int unchanged_borders = 0, curr_border, num_borders int borders[6] @@ -649,7 +717,6 @@ def _compute_thin_image(pixel_type[:, :, ::1] img not None): bint no_change list simple_border_points pixel_type[::1] neighb = np.zeros(27, dtype=np.uint8) - borders[:] = [4, 3, 2, 1, 5, 6] # no need to worry about the z direction if the original image is 2D. diff --git a/skimage/morphology/skel.py b/skimage/morphology/skel.py index 2da3979e..5af021b9 100644 --- a/skimage/morphology/skel.py +++ b/skimage/morphology/skel.py @@ -27,6 +27,7 @@ def _prepare_image(img_in): def _postprocess_image(img_o): """Clip the image (padding is an implementation detail), convert to b/w. + If the original was 2D, convert back to 2D. """ img_oo = img_o[1:-1, 1:-1, 1:-1] img_oo = img_oo.squeeze() @@ -35,6 +36,8 @@ def _postprocess_image(img_o): def compute_thin_image(img_in): + """Compute the thin image. + """ img = _prepare_image(img_in) img = np.asarray(_compute_thin_image(img)) img = _postprocess_image(img)