From 06d5a91f57dc7ab4c92c251e253a64ee018bdabf Mon Sep 17 00:00:00 2001 From: Evgeni Burovski Date: Tue, 26 Jan 2016 19:48:08 +0000 Subject: [PATCH] MAINT: skel3d: move computations to cython --- skimage/morphology/_skel.pyx | 638 +++++++++++++++++++++++++++++++++++ skimage/morphology/setup.py | 3 + skimage/morphology/skel.py | 637 +--------------------------------- 3 files changed, 644 insertions(+), 634 deletions(-) create mode 100644 skimage/morphology/_skel.pyx diff --git a/skimage/morphology/_skel.pyx b/skimage/morphology/_skel.pyx new file mode 100644 index 00000000..3c24bca1 --- /dev/null +++ b/skimage/morphology/_skel.pyx @@ -0,0 +1,638 @@ +from __future__ import division, print_function, absolute_import + +import numpy as np + + +def get_neighborhood(img, p, r, c): + """Get the neighborhood of a pixel. + + Assume zero boundary conditions. Image is already padded, so no + out-of-bounds checking. + """ + neighborhood = np.zeros(27, dtype=np.uint8) + + 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] + + return neighborhood + + +###### look-up tables +def fill_numpoints_LUT(n=256): + p = int(np.log2(n) + 1) + return np.sum(np.arange(n)[:, None] & (1 << np.arange(p)) != 0, axis=1) + +NUMPOINTS_LUT = fill_numpoints_LUT() + + +def fill_Euler_LUT(): + LUT = np.zeros(256, dtype=int) + + LUT[1] = 1 + LUT[3] = -1 + LUT[5] = -1 + LUT[7] = 1 + LUT[9] = -3 + LUT[11] = -1 + LUT[13] = -1 + LUT[15] = 1 + LUT[17] = -1 + LUT[19] = 1 + LUT[21] = 1 + LUT[23] = -1 + LUT[25] = 3 + LUT[27] = 1 + LUT[29] = 1 + LUT[31] = -1 + LUT[33] = -3 + LUT[35] = -1 + LUT[37] = 3 + LUT[39] = 1 + LUT[41] = 1 + LUT[43] = -1 + LUT[45] = 3 + LUT[47] = 1 + LUT[49] = -1 + LUT[51] = 1 + + LUT[53] = 1 + LUT[55] = -1 + LUT[57] = 3 + LUT[59] = 1 + LUT[61] = 1 + LUT[63] = -1 + LUT[65] = -3 + LUT[67] = 3 + LUT[69] = -1 + LUT[71] = 1 + LUT[73] = 1 + LUT[75] = 3 + LUT[77] = -1 + LUT[79] = 1 + LUT[81] = -1 + LUT[83] = 1 + LUT[85] = 1 + LUT[87] = -1 + LUT[89] = 3 + LUT[91] = 1 + LUT[93] = 1 + LUT[95] = -1 + LUT[97] = 1 + LUT[99] = 3 + LUT[101] = 3 + LUT[103] = 1 + + LUT[105] = 5 + LUT[107] = 3 + LUT[109] = 3 + LUT[111] = 1 + LUT[113] = -1 + LUT[115] = 1 + LUT[117] = 1 + LUT[119] = -1 + LUT[121] = 3 + LUT[123] = 1 + LUT[125] = 1 + LUT[127] = -1 + LUT[129] = -7 + LUT[131] = -1 + LUT[133] = -1 + LUT[135] = 1 + LUT[137] = -3 + LUT[139] = -1 + LUT[141] = -1 + LUT[143] = 1 + LUT[145] = -1 + LUT[147] = 1 + LUT[149] = 1 + LUT[151] = -1 + LUT[153] = 3 + LUT[155] = 1 + + LUT[157] = 1 + LUT[159] = -1 + LUT[161] = -3 + LUT[163] = -1 + LUT[165] = 3 + LUT[167] = 1 + LUT[169] = 1 + LUT[171] = -1 + LUT[173] = 3 + LUT[175] = 1 + LUT[177] = -1 + LUT[179] = 1 + LUT[181] = 1 + LUT[183] = -1 + LUT[185] = 3 + LUT[187] = 1 + LUT[189] = 1 + LUT[191] = -1 + LUT[193] = -3 + LUT[195] = 3 + LUT[197] = -1 + LUT[199] = 1 + LUT[201] = 1 + LUT[203] = 3 + LUT[205] = -1 + LUT[207] = 1 + + LUT[209] = -1 + LUT[211] = 1 + LUT[213] = 1 + LUT[215] = -1 + LUT[217] = 3 + LUT[219] = 1 + LUT[221] = 1 + LUT[223] = -1 + LUT[225] = 1 + LUT[227] = 3 + LUT[229] = 3 + LUT[231] = 1 + LUT[233] = 5 + LUT[235] = 3 + LUT[237] = 3 + LUT[239] = 1 + LUT[241] = -1 + LUT[243] = 1 + LUT[245] = 1 + LUT[247] = -1 + LUT[249] = 3 + LUT[251] = 1 + LUT[253] = 1 + LUT[255] = -1 + return LUT + +LUT = fill_Euler_LUT() + + +### Octants (indexOctantXXX functions) +OCTANTS = tuple(range(8)) +NEB, NWB, SEB, SWB, NEU, NWU, SEU, SWU = OCTANTS + +neib_idx = np.empty((8, 7), dtype=int) +neib_idx[NEB, ...] = [2, 1, 11, 10, 5, 4, 14] +neib_idx[NWB, ...] = [0, 9, 3, 12, 1, 10, 4] +neib_idx[SEB, ...] = [8, 7, 17, 16, 5, 4, 14] +neib_idx[SWB, ...] = [6, 15, 7, 16, 3, 12, 4] +neib_idx[NEU, ...] = [20, 23, 19, 22, 11, 14, 10] +neib_idx[NWU, ...] = [18, 21, 9, 12, 19, 22, 10] +neib_idx[SEU, ...] = [26, 23, 17, 14, 25, 22, 16] +neib_idx[SWU, ...] = [24, 25, 15, 16, 21, 22, 12] + +def index_octants(octant, neighbors): + n = 1 + for j, idx in enumerate(neib_idx[octant]): + if neighbors[idx] == 1: + n |= 2**(7 - j) + return n + + +def is_surfacepoint(neighbors, points_LUT): + for octant in OCTANTS: + n = index_octants(octant, neighbors) + if n not in (240, 165, 170) and points_LUT[n] > 2: + return False + return True + + +def is_Euler_invariant(neighbors): + """Check if a point is Euler invariant. + + Calculate Euler characteristc for each octant and sum up. + + Parameters + ---------- + neighbors : ndarray, shape (27,) + neighbors of a point + + Returns + ------- + bool + + """ + euler_char = 0 + for octant in OCTANTS: + n = index_octants(octant, neighbors) + euler_char += LUT[n] + return euler_char == 0 + + +def is_simple_point(neighbors): + """Check is a point is a Simple Point. + + 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. + + Parameters + ---------- + neighbors : ndarray, 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]) + cube = np.r_[neighbors[:13], neighbors[14:]] + + # set initial label + label = 2 + + # 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) + else: + raise ValueError("Never be here. i = %s" % i) + label += 1 + if label - 2 >= 2: + return False + return True + + +def octree_labeling(octant, label, 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. + + Parameters + ---------- + octant : int + octant index + label : int + the current label of the center point + cube : ndarray, shape(26,) + local neighborhood of the point + + """ + # check if there are points in the octant with value 1 + if octant == 1: + # set points in this octant to current label + # and recursive labeling of adjacent octants + if cube[0] == 1: + cube[0] = label + if cube[1] == 1: + cube[1] = label + octree_labeling(2, label, cube) + if cube[3] == 1: + cube[3] = label + octree_labeling(3, label, cube) + if cube[4] == 1: + cube[4] = label + octree_labeling(2, label, cube) + octree_labeling(3, label, cube) + octree_labeling(4, label, cube) + if cube[9] == 1: + cube[9] = label + octree_labeling(5, label, cube) + if cube[10] == 1: + cube[10] = label + octree_labeling(2, label, cube) + octree_labeling(5, label, cube) + octree_labeling(6, label, cube) + if cube[12] == 1: + cube[12] = label + octree_labeling(3, label, cube) + octree_labeling(5, label, cube) + octree_labeling(7, label, cube) + + if octant == 2: + if cube[1] == 1: + cube[1] = label + octree_labeling(1, label, cube) + if cube[4] == 1: + cube[4] = label + octree_labeling(1, label, cube) + octree_labeling(3, label, cube) + octree_labeling(4, label, cube) + if cube[10] == 1: + cube[10] = label + octree_labeling(1, label, cube) + octree_labeling(5, label, cube) + octree_labeling(6, label, cube) + if cube[2] == 1: + cube[2] = label + if cube[5] == 1: + cube[5] = label + octree_labeling(4, label, cube) + if cube[11] == 1: + cube[11] = label + octree_labeling(6, label, cube) + if cube[13] == 1: + cube[13] = label + octree_labeling(4, label, cube) + octree_labeling(6, label, cube) + octree_labeling(8, label, cube) + + if octant ==3: + if cube[3] == 1: + cube[3] = label + octree_labeling(1, label, cube) + if cube[4] == 1: + cube[4] = label + octree_labeling(1, label, cube) + octree_labeling(2, label, cube) + octree_labeling(4, label, cube) + if cube[12] == 1: + cube[12] = label + octree_labeling(1, label, cube) + octree_labeling(5, label, cube) + octree_labeling(7, label, cube) + if cube[6] == 1: + cube[6] = label + if cube[7] == 1: + cube[7] = label + octree_labeling(4, label, cube) + if cube[14] == 1: + cube[14] = label + octree_labeling(7, label, cube) + if cube[15] == 1: + cube[15] = label + octree_labeling(4, label, cube) + octree_labeling(7, label, cube) + octree_labeling(8, label, cube) + + if octant == 4: + if cube[4] == 1: + cube[4] = label + octree_labeling(1, label, cube) + octree_labeling(2, label, cube) + octree_labeling(3, label, cube) + if cube[5] == 1: + cube[5] = label + octree_labeling(2, label, cube) + if cube[13] == 1: + cube[13] = label + octree_labeling(2, label, cube) + octree_labeling(6, label, cube) + octree_labeling(8, label, cube) + if cube[7] == 1: + cube[7] = label + octree_labeling(3, label, cube) + if cube[15] == 1: + cube[15] = label + octree_labeling(3, label, cube) + octree_labeling(7, label, cube) + octree_labeling(8, label, cube) + if cube[8] == 1: + cube[8] = label + if cube[16] == 1: + cube[16] = label + octree_labeling(8, label, cube) + + if octant == 5: + if cube[9] == 1: + cube[9] = label + octree_labeling(1, label, cube) + if cube[10] == 1: + cube[10] = label + octree_labeling(1, label, cube) + octree_labeling(2, label, cube) + octree_labeling(6, label, cube) + if cube[12] == 1: + cube[12] = label + octree_labeling(1, label, cube) + octree_labeling(3, label, cube) + octree_labeling(7, label, cube) + if cube[17] == 1: + cube[17] = label + if cube[18] == 1: + cube[18] = label + octree_labeling(6, label, cube) + if cube[20] == 1: + cube[20] = label + octree_labeling(7, label, cube) + if cube[21] == 1: + cube[21] = label + octree_labeling(6, label, cube) + octree_labeling(7, label, cube) + octree_labeling(8, label, cube) + + if octant == 6: + if cube[10] == 1: + cube[10] = label + octree_labeling(1, label, cube) + octree_labeling(2, label, cube) + octree_labeling(5, label, cube) + if cube[11] == 1: + cube[11] = label + octree_labeling(2, label, cube) + if cube[13] == 1: + cube[13] = label + octree_labeling(2, label, cube) + octree_labeling(4, label, cube) + octree_labeling(8, label, cube) + if cube[18] == 1: + cube[18] = label + octree_labeling(5, label, cube) + if cube[21] == 1: + cube[21] = label + octree_labeling(5, label, cube) + octree_labeling(7, label, cube) + octree_labeling(8, label, cube) + if cube[19] == 1: + cube[19] = label + if cube[22] == 1: + cube[22] = label + octree_labeling(8, label, cube) + + if octant == 7: + if cube[12] == 1: + cube[12] = label + octree_labeling(1, label, cube) + octree_labeling(3, label, cube) + octree_labeling(5, label, cube) + if cube[14] == 1: + cube[14] = label + octree_labeling(3, label, cube) + if cube[15] == 1: + cube[15] = label + octree_labeling(3, label, cube) + octree_labeling(4, label, cube) + octree_labeling(8, label, cube) + if cube[20] == 1: + cube[20] = label + octree_labeling(5, label, cube) + if cube[21] == 1: + cube[21] = label + octree_labeling(5, label, cube) + octree_labeling(6, label, cube) + octree_labeling(8, label, cube) + if cube[23] == 1: + cube[23] = label + if cube[24] == 1: + cube[24] = label + octree_labeling(8, label, cube) + + if octant == 8: + if cube[13] == 1: + cube[13] = label + octree_labeling(2, label, cube) + octree_labeling(4, label, cube) + octree_labeling(6, label, cube) + if cube[15] == 1: + cube[15] = label + octree_labeling(3, label, cube) + octree_labeling(4, label, cube) + octree_labeling(7, label, cube) + if cube[16] == 1: + cube[16] = label + octree_labeling(4, label, cube) + if cube[21] == 1: + cube[21] = label + octree_labeling(5, label, cube) + octree_labeling(6, label, cube) + octree_labeling(7, label, cube) + if cube[22] == 1: + cube[22] = label + octree_labeling(6, label, cube) + if cube[24] == 1: + cube[24] = label + octree_labeling(7, label, cube) + if cube[25] == 1: + cube[25] = label + + +def _loop_through(img, curr_border): + """Inner loop of compute_thin_image. + + return simple_border_points as a list to be rechecked sequentially. + """ + # loop through the image + # NB: each loop is from 1 to size-1: img is padded from all sides + simple_border_points = [] + + ### XXX: 2D images + ### if the original is 2D, img.shape[0] == 3, the algorithm removes too much + ### because all points are considered 'boundary' in the 3rd direction. + ### Hence just bail out + if img.shape[0] == 3 and curr_border in (5, 6): + print("skipping curr_border = ", curr_border) + return [] + + 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 + + neighborhood = get_neighborhood(img, p, r, c) + + # check if (p, r, c) is an endpoint. An endpoint has exactly + # one neighbor in the 26-neighborhood. + # The center pixel is counted, thus r.h.s. is 2 + if neighborhood.sum() == 2: + continue + + # check if point is Euler invariant (condition 1 in [Lee94]) + # if it is not, it's not deletable + if not is_Euler_invariant(neighborhood): + continue + + # check if point is simple (i.e., deletion does not + # change connectivity in the 3x3x3 neighborhood) + # this are conditions 2 and 3 in [Lee94] + if not is_simple_point(neighborhood): + continue + + # ok, add (p, r, c) to the list of simple border points + simple_border_points.append((p, r, c)) + return simple_border_points + + +def _compute_thin_image(img): + ### compute + unchanged_borders = 0 + + # loop through the image several times until there is no change for all + # the six border types + while unchanged_borders < 6: + unchanged_borders = 0 + for curr_border in (4, 3, 2, 1, 5, 6): + + simple_border_points = _loop_through(img, curr_border) + print(curr_border, " : ", simple_border_points, '\n') + + # sequential re-checking to preserve connectivity when deleting + # in a parallel way + no_change = True + for pt in simple_border_points: + p, r, c = pt + neighb = get_neighborhood(img, p, r, c) + if is_simple_point(neighb): + img[p, r, c] = 0 + no_change = False + else: + print(" *** ", pt, is_simple_point(neighb)) + + if no_change: + unchanged_borders += 1 + simple_border_points = [] + + return img diff --git a/skimage/morphology/setup.py b/skimage/morphology/setup.py index dbbcad8b..5a26ac52 100644 --- a/skimage/morphology/setup.py +++ b/skimage/morphology/setup.py @@ -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(['_skel.pyx'], 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('_skel', sources=['_skel.c'], + include_dirs=[get_numpy_include_dirs()]) return config diff --git a/skimage/morphology/skel.py b/skimage/morphology/skel.py index e53c7ee6..5f8d9fad 100644 --- a/skimage/morphology/skel.py +++ b/skimage/morphology/skel.py @@ -2,6 +2,8 @@ from __future__ import division, print_function, absolute_import import numpy as np +from ._skel import _compute_thin_image + def _prepare_image(img_in): """Convert to a binary image, pad the it w/ zeros, and ensure it's 3D. @@ -34,642 +36,9 @@ def _postprocess_image(img_o): return img_oo -def get_neighborhood(img, p, r, c): - """Get the neighborhood of a pixel. - - Assume zero boundary conditions. Image is already padded, so no - out-of-bounds checking. - """ - neighborhood = np.zeros(27, dtype=np.uint8) - - 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] - - return neighborhood - - -###### look-up tables -def fill_numpoints_LUT(n=256): - p = int(np.log2(n) + 1) - return np.sum(np.arange(n)[:, None] & (1 << np.arange(p)) != 0, axis=1) - -NUMPOINTS_LUT = fill_numpoints_LUT() - - -def fill_Euler_LUT(): - LUT = np.zeros(256, dtype=int) - - LUT[1] = 1 - LUT[3] = -1 - LUT[5] = -1 - LUT[7] = 1 - LUT[9] = -3 - LUT[11] = -1 - LUT[13] = -1 - LUT[15] = 1 - LUT[17] = -1 - LUT[19] = 1 - LUT[21] = 1 - LUT[23] = -1 - LUT[25] = 3 - LUT[27] = 1 - LUT[29] = 1 - LUT[31] = -1 - LUT[33] = -3 - LUT[35] = -1 - LUT[37] = 3 - LUT[39] = 1 - LUT[41] = 1 - LUT[43] = -1 - LUT[45] = 3 - LUT[47] = 1 - LUT[49] = -1 - LUT[51] = 1 - - LUT[53] = 1 - LUT[55] = -1 - LUT[57] = 3 - LUT[59] = 1 - LUT[61] = 1 - LUT[63] = -1 - LUT[65] = -3 - LUT[67] = 3 - LUT[69] = -1 - LUT[71] = 1 - LUT[73] = 1 - LUT[75] = 3 - LUT[77] = -1 - LUT[79] = 1 - LUT[81] = -1 - LUT[83] = 1 - LUT[85] = 1 - LUT[87] = -1 - LUT[89] = 3 - LUT[91] = 1 - LUT[93] = 1 - LUT[95] = -1 - LUT[97] = 1 - LUT[99] = 3 - LUT[101] = 3 - LUT[103] = 1 - - LUT[105] = 5 - LUT[107] = 3 - LUT[109] = 3 - LUT[111] = 1 - LUT[113] = -1 - LUT[115] = 1 - LUT[117] = 1 - LUT[119] = -1 - LUT[121] = 3 - LUT[123] = 1 - LUT[125] = 1 - LUT[127] = -1 - LUT[129] = -7 - LUT[131] = -1 - LUT[133] = -1 - LUT[135] = 1 - LUT[137] = -3 - LUT[139] = -1 - LUT[141] = -1 - LUT[143] = 1 - LUT[145] = -1 - LUT[147] = 1 - LUT[149] = 1 - LUT[151] = -1 - LUT[153] = 3 - LUT[155] = 1 - - LUT[157] = 1 - LUT[159] = -1 - LUT[161] = -3 - LUT[163] = -1 - LUT[165] = 3 - LUT[167] = 1 - LUT[169] = 1 - LUT[171] = -1 - LUT[173] = 3 - LUT[175] = 1 - LUT[177] = -1 - LUT[179] = 1 - LUT[181] = 1 - LUT[183] = -1 - LUT[185] = 3 - LUT[187] = 1 - LUT[189] = 1 - LUT[191] = -1 - LUT[193] = -3 - LUT[195] = 3 - LUT[197] = -1 - LUT[199] = 1 - LUT[201] = 1 - LUT[203] = 3 - LUT[205] = -1 - LUT[207] = 1 - - LUT[209] = -1 - LUT[211] = 1 - LUT[213] = 1 - LUT[215] = -1 - LUT[217] = 3 - LUT[219] = 1 - LUT[221] = 1 - LUT[223] = -1 - LUT[225] = 1 - LUT[227] = 3 - LUT[229] = 3 - LUT[231] = 1 - LUT[233] = 5 - LUT[235] = 3 - LUT[237] = 3 - LUT[239] = 1 - LUT[241] = -1 - LUT[243] = 1 - LUT[245] = 1 - LUT[247] = -1 - LUT[249] = 3 - LUT[251] = 1 - LUT[253] = 1 - LUT[255] = -1 - return LUT - -LUT = fill_Euler_LUT() - - -### Octants (indexOctantXXX functions) -OCTANTS = tuple(range(8)) -NEB, NWB, SEB, SWB, NEU, NWU, SEU, SWU = OCTANTS - -neib_idx = np.empty((8, 7), dtype=int) -neib_idx[NEB, ...] = [2, 1, 11, 10, 5, 4, 14] -neib_idx[NWB, ...] = [0, 9, 3, 12, 1, 10, 4] -neib_idx[SEB, ...] = [8, 7, 17, 16, 5, 4, 14] -neib_idx[SWB, ...] = [6, 15, 7, 16, 3, 12, 4] -neib_idx[NEU, ...] = [20, 23, 19, 22, 11, 14, 10] -neib_idx[NWU, ...] = [18, 21, 9, 12, 19, 22, 10] -neib_idx[SEU, ...] = [26, 23, 17, 14, 25, 22, 16] -neib_idx[SWU, ...] = [24, 25, 15, 16, 21, 22, 12] - -def index_octants(octant, neighbors): - n = 1 - for j, idx in enumerate(neib_idx[octant]): - if neighbors[idx] == 1: - n |= 2**(7 - j) - return n - - -def is_surfacepoint(neighbors, points_LUT): - for octant in OCTANTS: - n = index_octants(octabt, neighbors) - if n not in (240, 165, 170) and points_LUT[n] > 2: - return False - return True - - -def is_Euler_invariant(neighbors): - """Check if a point is Euler invariant. - - Calculate Euler characteristc for each octant and sum up. - - Parameters - ---------- - neighbors : ndarray, shape (27,) - neighbors of a point - - Returns - ------- - bool - - """ - euler_char = 0 - for octant in OCTANTS: - n = index_octants(octant, neighbors) - euler_char += LUT[n] - return euler_char == 0 - - -def is_simple_point(neighbors): - """Check is a point is a Simple Point. - - 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. - - Parameters - ---------- - neighbors : ndarray, 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]) - cube = np.r_[neighbors[:13], neighbors[14:]] - - # set initial label - label = 2 - - # 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) - else: - raise ValueError("Never be here. i = %s" % i) - label += 1 - if label - 2 >= 2: - return False - return True - - -def octree_labeling(octant, label, 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. - - Parameters - ---------- - octant : int - octant index - label : int - the current label of the center point - cube : ndarray, shape(26,) - local neighborhood of the point - - """ - # check if there are points in the octant with value 1 - if octant == 1: - # set points in this octant to current label - # and recursive labeling of adjacent octants - if cube[0] == 1: - cube[0] = label - if cube[1] == 1: - cube[1] = label - octree_labeling(2, label, cube) - if cube[3] == 1: - cube[3] = label - octree_labeling(3, label, cube) - if cube[4] == 1: - cube[4] = label - octree_labeling(2, label, cube) - octree_labeling(3, label, cube) - octree_labeling(4, label, cube) - if cube[9] == 1: - cube[9] = label - octree_labeling(5, label, cube) - if cube[10] == 1: - cube[10] = label - octree_labeling(2, label, cube) - octree_labeling(5, label, cube) - octree_labeling(6, label, cube) - if cube[12] == 1: - cube[12] = label - octree_labeling(3, label, cube) - octree_labeling(5, label, cube) - octree_labeling(7, label, cube) - - if octant == 2: - if cube[1] == 1: - cube[1] = label - octree_labeling(1, label, cube) - if cube[4] == 1: - cube[4] = label - octree_labeling(1, label, cube) - octree_labeling(3, label, cube) - octree_labeling(4, label, cube) - if cube[10] == 1: - cube[10] = label - octree_labeling(1, label, cube) - octree_labeling(5, label, cube) - octree_labeling(6, label, cube) - if cube[2] == 1: - cube[2] = label - if cube[5] == 1: - cube[5] = label - octree_labeling(4, label, cube) - if cube[11] == 1: - cube[11] = label - octree_labeling(6, label, cube) - if cube[13] == 1: - cube[13] = label - octree_labeling(4, label, cube) - octree_labeling(6, label, cube) - octree_labeling(8, label, cube) - - if octant ==3: - if cube[3] == 1: - cube[3] = label - octree_labeling(1, label, cube) - if cube[4] == 1: - cube[4] = label - octree_labeling(1, label, cube) - octree_labeling(2, label, cube) - octree_labeling(4, label, cube) - if cube[12] == 1: - cube[12] = label - octree_labeling(1, label, cube) - octree_labeling(5, label, cube) - octree_labeling(7, label, cube) - if cube[6] == 1: - cube[6] = label - if cube[7] == 1: - cube[7] = label - octree_labeling(4, label, cube) - if cube[14] == 1: - cube[14] = label - octree_labeling(7, label, cube) - if cube[15] == 1: - cube[15] = label - octree_labeling(4, label, cube) - octree_labeling(7, label, cube) - octree_labeling(8, label, cube) - - if octant == 4: - if cube[4] == 1: - cube[4] = label - octree_labeling(1, label, cube) - octree_labeling(2, label, cube) - octree_labeling(3, label, cube) - if cube[5] == 1: - cube[5] = label - octree_labeling(2, label, cube) - if cube[13] == 1: - cube[13] = label - octree_labeling(2, label, cube) - octree_labeling(6, label, cube) - octree_labeling(8, label, cube) - if cube[7] == 1: - cube[7] = label - octree_labeling(3, label, cube) - if cube[15] == 1: - cube[15] = label - octree_labeling(3, label, cube) - octree_labeling(7, label, cube) - octree_labeling(8, label, cube) - if cube[8] == 1: - cube[8] = label - if cube[16] == 1: - cube[16] = label - octree_labeling(8, label, cube) - - if octant == 5: - if cube[9] == 1: - cube[9] = label - octree_labeling(1, label, cube) - if cube[10] == 1: - cube[10] = label - octree_labeling(1, label, cube) - octree_labeling(2, label, cube) - octree_labeling(6, label, cube) - if cube[12] == 1: - cube[12] = label - octree_labeling(1, label, cube) - octree_labeling(3, label, cube) - octree_labeling(7, label, cube) - if cube[17] == 1: - cube[17] = label - if cube[18] == 1: - cube[18] = label - octree_labeling(6, label, cube) - if cube[20] == 1: - cube[20] = label - octree_labeling(7, label, cube) - if cube[21] == 1: - cube[21] = label - octree_labeling(6, label, cube) - octree_labeling(7, label, cube) - octree_labeling(8, label, cube) - - if octant == 6: - if cube[10] == 1: - cube[10] = label - octree_labeling(1, label, cube) - octree_labeling(2, label, cube) - octree_labeling(5, label, cube) - if cube[11] == 1: - cube[11] = label - octree_labeling(2, label, cube) - if cube[13] == 1: - cube[13] = label - octree_labeling(2, label, cube) - octree_labeling(4, label, cube) - octree_labeling(8, label, cube) - if cube[18] == 1: - cube[18] = label - octree_labeling(5, label, cube) - if cube[21] == 1: - cube[21] = label - octree_labeling(5, label, cube) - octree_labeling(7, label, cube) - octree_labeling(8, label, cube) - if cube[19] == 1: - cube[19] = label - if cube[22] == 1: - cube[22] = label - octree_labeling(8, label, cube) - - if octant == 7: - if cube[12] == 1: - cube[12] = label - octree_labeling(1, label, cube) - octree_labeling(3, label, cube) - octree_labeling(5, label, cube) - if cube[14] == 1: - cube[14] = label - octree_labeling(3, label, cube) - if cube[15] == 1: - cube[15] = label - octree_labeling(3, label, cube) - octree_labeling(4, label, cube) - octree_labeling(8, label, cube) - if cube[20] == 1: - cube[20] = label - octree_labeling(5, label, cube) - if cube[21] == 1: - cube[21] = label - octree_labeling(5, label, cube) - octree_labeling(6, label, cube) - octree_labeling(8, label, cube) - if cube[23] == 1: - cube[23] = label - if cube[24] == 1: - cube[24] = label - octree_labeling(8, label, cube) - - if octant == 8: - if cube[13] == 1: - cube[13] = label - octree_labeling(2, label, cube) - octree_labeling(4, label, cube) - octree_labeling(6, label, cube) - if cube[15] == 1: - cube[15] = label - octree_labeling(3, label, cube) - octree_labeling(4, label, cube) - octree_labeling(7, label, cube) - if cube[16] == 1: - cube[16] = label - octree_labeling(4, label, cube) - if cube[21] == 1: - cube[21] = label - octree_labeling(5, label, cube) - octree_labeling(6, label, cube) - octree_labeling(7, label, cube) - if cube[22] == 1: - cube[22] = label - octree_labeling(6, label, cube) - if cube[24] == 1: - cube[24] = label - octree_labeling(7, label, cube) - if cube[25] == 1: - cube[25] = label - - -def _loop_through(img, curr_border): - """Inner loop of compute_thin_image. - - return simple_border_points as a list to be rechecked sequentially. - """ - # loop through the image - # NB: each loop is from 1 to size-1: img is padded from all sides - simple_border_points = [] - - ### XXX: 2D images - ### if the original is 2D, img.shape[0] == 3, the algorithm removes too much - ### because all points are considered 'boundary' in the 3rd direction. - ### Hence just bail out - if img.shape[0] == 3 and curr_border in (5, 6): - print("skipping curr_border = ", curr_border) - return [] - - 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 - - neighborhood = get_neighborhood(img, p, r, c) - - # check if (p, r, c) is an endpoint. An endpoint has exactly - # one neighbor in the 26-neighborhood. - # The center pixel is counted, thus r.h.s. is 2 - if neighborhood.sum() == 2: - continue - - # check if point is Euler invariant (condition 1 in [Lee94]) - # if it is not, it's not deletable - if not is_Euler_invariant(neighborhood): - continue - - # check if point is simple (i.e., deletion does not - # change connectivity in the 3x3x3 neighborhood) - # this are conditions 2 and 3 in [Lee94] - if not is_simple_point(neighborhood): - continue - - # ok, add (p, r, c) to the list of simple border points - simple_border_points.append((p, r, c)) - return simple_border_points - - def compute_thin_image(img_in): - - ### prepare img = _prepare_image(img_in) - - ### compute - unchanged_borders = 0 - - # loop through the image several times until there is no change for all - # the six border types - while unchanged_borders < 6: - unchanged_borders = 0 - for curr_border in (4, 3, 2, 1, 5, 6): - - simple_border_points = _loop_through(img, curr_border) - print(curr_border, " : ", simple_border_points, '\n') - - # sequential re-checking to preserve connectivity when deleting - # in a parallel way - no_change = True - for pt in simple_border_points: - p, r, c = pt - neighb = get_neighborhood(img, p, r, c) - if is_simple_point(neighb): - img[p, r, c] = 0 - no_change = False - else: - print(" *** ", pt, is_simple_point(neighb)) - - if no_change: - unchanged_borders += 1 - simple_border_points = [] - + img = _compute_thin_image(img) img = _postprocess_image(img) return img