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

MAINT: skel3d: address review comments

Browse Source 
* simplify filling in the neighb_idx LUT
* underscore references
* preserve the color range of the original image
This commit is contained in:
Evgeni Burovski
2016-02-20 17:26:12 +00:00
parent 456a5252c9
commit 85f46c6fe8
3 changed files with 32 additions and 111 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/morphology/_skel.pyx
+25 -33
View File
@@ -1,8 +1,8 @@
"""
This is an implementation of the 2D/3D thinning algorithm
of [Lee94] of binary images, based on [IAC15].
of [Lee94]_ of binary images, based on [IAC15]_.
The original Java code [IAC15] carries the following message:
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
@@ -49,7 +49,7 @@ def _compute_thin_image(pixel_type[:, :, ::1] img not None):
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,
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.
@@ -112,9 +112,9 @@ cdef list _loop_through(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
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.
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.
@@ -160,14 +160,14 @@ cdef list _loop_through(pixel_type[:, :, ::1] img,
if is_endpoint(neighborhood):
continue
# check if point is Euler invariant (condition 1 in [Lee94]):
# check if point is Euler invariant (condition 1 in [Lee94]_):
# if it is not, it's not deletable.
if not is_Euler_invariant(neighborhood, Euler_LUT, neighb_idx):
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]
# this are conditions 2 and 3 in [Lee94]_
if not is_simple_point(neighborhood):
continue
@@ -186,9 +186,9 @@ cdef void get_neighborhood(pixel_type[:, :, ::1] img,
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
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].
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]
@@ -231,7 +231,7 @@ cdef void get_neighborhood(pixel_type[:, :, ::1] img,
def fill_Euler_LUT():
""" Look-up table for preserving Euler characteristic.
This is column $\delta G_{26}$ of Table 2 of [Lee94].
This is column $\delta G_{26}$ of Table 2 of [Lee94]_.
"""
LUT = np.zeros(256, dtype=np.intc)
@@ -371,25 +371,17 @@ def fill_Euler_LUT():
cdef int[::1] LUT = fill_Euler_LUT()
### Octants (indexOctantXXX functions)
def fill_neighbor_idx():
"""Fill the look-up table for indexing octants for computing the Euler
characteristics.
See index_octants and is_Euler_invariant routines below.
"""
NEB, NWB, SEB, SWB, NEU, NWU, SEU, SWU = tuple(range(8))
_neighb_idx = np.empty((8, 7), dtype=np.intc)
_neighb_idx[NEB, ...] = [2, 1, 11, 10, 5, 4, 14]
_neighb_idx[NWB, ...] = [0, 9, 3, 12, 1, 10, 4]
_neighb_idx[SEB, ...] = [8, 7, 17, 16, 5, 4, 14]
_neighb_idx[SWB, ...] = [6, 15, 7, 16, 3, 12, 4]
_neighb_idx[NEU, ...] = [20, 23, 19, 22, 11, 14, 10]
_neighb_idx[NWU, ...] = [18, 21, 9, 12, 19, 22, 10]
_neighb_idx[SEU, ...] = [26, 23, 17, 14, 25, 22, 16]
_neighb_idx[SWU, ...] = [24, 25, 15, 16, 21, 22, 12]
return _neighb_idx
cdef int[:, ::1] NEIGHB_IDX = fill_neighbor_idx()
# Fill the look-up table for indexing octants for computing the Euler
# characteristic. See index_octants and is_Euler_invariant routines below.
cdef int[:, ::1] NEIGHB_IDX = np.array([[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
], dtype=np.intc)
@cython.boundscheck(False)
@@ -402,7 +394,7 @@ cdef int index_octants(int octant,
for j in range(7):
idx = neib_idx[octant, j]
if neighbors[idx] == 1:
n |= 2 ** (7 - j) # XXX hardcode powers?
n |= 1 << (7 - j) # XXX hardcode powers?
return n
@@ -449,7 +441,7 @@ cdef bint is_Euler_invariant(pixel_type neighbors[],
cdef bint is_simple_point(pixel_type 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.
@@ -465,7 +457,7 @@ cdef bint is_simple_point(pixel_type neighbors[]):
"""
# copy neighbors for labeling
# ignore center pixel (i=13) when counting (see [Lee94])
# 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))
@@ -509,7 +501,7 @@ cdef void octree_labeling(int octant, int label, pixel_type 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.
See Figs. 6 and 7 of [Lee94]_ for the values of indices.
Parameters
----------
skimage/morphology/skel.py
+7 -5
View File
@@ -15,6 +15,8 @@ def _prepare_image(img_in):
if img.ndim == 2:
img = img.reshape((1,) + img.shape)
intensity = img.max()
# normalize to binary
img[img != 0] = 1
@@ -22,25 +24,25 @@ def _prepare_image(img_in):
img_o = np.zeros(tuple(s + 2 for s in img.shape),
dtype=np.uint8)
img_o[1:-1, 1:-1, 1:-1] = img.astype(np.uint8)
return img_o
return img_o, intensity
def _postprocess_image(img_o):
def _postprocess_image(img_o, intensity):
"""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()
img_oo *= 255
img_oo *= intensity
return img_oo
def compute_thin_image(img_in):
"""Compute the thin image.
"""
img = _prepare_image(img_in)
img, intensity = _prepare_image(img_in)
img = np.asarray(_compute_thin_image(img))
img = _postprocess_image(img)
img = _postprocess_image(img, intensity)
return img
skimage/morphology/trivial.py
-73
View File
@@ -1,73 +0,0 @@
import numpy as np
import matplotlib.pyplot as plt
from skel import prepare_image, compute_thin_image
def get_rhombus(n=64, L=22, width=7):
img = np.zeros((n, n), dtype=np.int8)
x = np.arange(L, dtype=int)
y = L - x
for w in range(width):
img[x + n//2, y + n//2 + w] = 1
img[-x + n//2, y + n//2 + w] = 1
img[x + n//2, -y + n//2 + w] = 1
img[-x + n//2, -y + n//2 + w] = 1
return img
def get_strip():
img = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=np.uint8)
return img
def get_square(n=64, L=16, width=4):
img = np.zeros((n, n), dtype=np.int8)
x = np.arange(-L//2 + 1, L//2)
for w in range(width):
img[n//2 + x, n//2 + L//2 + w] = 1
img[n//2 + x, n//2 - L//2 - w] = 1
img[n//2 + L//2 - w, n//2 + x] = 1
img[n//2 + x, n//2 - L//2 - w] = 1
return img
def get_loop():
img = np.loadtxt('Untitled.txt', dtype=np.uint8)
return img
if __name__ == "__main__":
## img = get_rhombus()
## img = get_strip()
img = get_loop()
x, y = np.nonzero(img)
plt.scatter(x, y, marker='s', color='b', s=40, alpha=0.3)
# skeletonize
img1 = prepare_image(img)
img1 = compute_thin_image(img1)
img1_2d = img1[1, 1:, 1:]
x, y = np.nonzero(img1_2d)
plt.scatter(x, y, marker='o', color='r')
plt.show()