scikit-image

1
#cython: cdivision=True
2
#cython: boundscheck=False
3
#cython: nonecheck=False
4
#cython: wraparound=False
5

6
import numpy as np
7
cimport numpy as cnp
8
cnp.import_array()
9

10

11
def _fast_skeletonize(cnp.uint8_t [:, ::1] image):
12
    """Optimized parts of the Zhang-Suen [1]_ skeletonization.
13
    Iteratively, pixels meeting removal criteria are removed,
14
    till only the skeleton remains (that is, no further removable pixel
15
    was found).
16

17
    Performs a hard-coded correlation to assign every neighborhood of 8 a
18
    unique number, which in turn is used in conjunction with a look up
19
    table to select the appropriate thinning criteria.
20

21
    Parameters
22
    ----------
23
    image : numpy.ndarray
24
        A binary image containing the objects to be skeletonized. '1'
25
        represents foreground, and '0' represents background.
26

27
    Returns
28
    -------
29
    skeleton : ndarray
30
        A matrix containing the thinned image.
31

32
    References
33
    ----------
34
    .. [1] A fast parallel algorithm for thinning digital patterns,
35
           T. Y. Zhang and C. Y. Suen, Communications of the ACM,
36
           March 1984, Volume 27, Number 3.
37

38
    """
39

40
    # look up table - there is one entry for each of the 2^8=256 possible
41
    # combinations of 8 binary neighbors. 1's, 2's and 3's are candidates
42
    # for removal at each iteration of the algorithm.
43
    cdef cnp.uint8_t *lut = [0, 0, 0, 1, 0, 0, 1, 3, 0, 0, 3, 1, 1, 0,
44
                             1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,
45
                             3, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0,
46
                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
47
                             2, 0, 0, 0, 3, 0, 2, 2, 0, 0, 0, 0, 0, 0,
48
                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
49
                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,
50
                             0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 2, 0, 0, 0,
51
                             3, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 3, 0,
52
                             2, 0, 0, 0, 3, 1, 0, 0, 1, 3, 0, 0, 0, 0,
53
                             0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
54
                             0, 0, 0, 0, 0, 1, 3, 1, 0, 0, 0, 0, 0, 0,
55
                             0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0,
56
                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 1, 3,
57
                             0, 0, 1, 3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
58
                             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
59
                             2, 3, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
60
                             0, 0, 3, 3, 0, 1, 0, 0, 0, 0, 2, 2, 0, 0,
61
                             2, 0, 0, 0]
62

63
    cdef cnp.uint8_t first_pass, pixel_removed, neighbors
64

65
    # indices for fast iteration
66
    cdef Py_ssize_t row, col, nrows, ncols, pass_num
67
    nrows, ncols = image.shape[:2]
68
    nrows += 2
69
    ncols += 2
70

71
    # Copy over the image into a larger version with a single pixel border
72
    # this removes the need to handle border cases below
73
    cdef cnp.uint8_t [:, ::1] skeleton = np.zeros((nrows, ncols),
74
                                                  dtype=np.uint8)
75
    skeleton[1:-1, 1:-1] = image
76
    cdef cnp.uint8_t [:, ::1] cleaned_skeleton = skeleton.copy()
77

78
    pixel_removed = True
79

80
    # the algorithm reiterates the thinning till
81
    # no further thinning occurred (variable pixel_removed set)
82
    with nogil:
83
        while pixel_removed:
84
            pixel_removed = False
85

86
            # there are two phases, in the first phase, pixels labeled
87
            # (see below) 1 and 3 are removed, in the second 2 and 3
88

89
            # nogil can't iterate through `(True, False)` because it is a Python
90
            # tuple. Use the fact that 0 is Falsy, and 1 is truthy in C
91
            # for the iteration instead.
92
            # for first_pass in (True, False):
93
            for pass_num in range(2):
94
                first_pass = (pass_num == 0)
95
                for row in range(1, nrows-1):
96
                    for col in range(1, ncols-1):
97
                        # all set pixels ...
98

99
                        if skeleton[row, col]:
100
                            # are correlated with a kernel
101
                            # (coefficients spread around here ...)
102
                            # to apply a unique number to every
103
                            # possible neighborhood ...
104

105
                            # which is used with the lut to find the
106
                            # "connectivity type"
107

108
                            neighbors = lut[skeleton[row - 1, col - 1] +
109
                                            2 * skeleton[row - 1, col] +
110
                                            4 * skeleton[row - 1, col + 1] +
111
                                            8 * skeleton[row, col + 1] +
112
                                            16 * skeleton[row + 1, col + 1] +
113
                                            32 * skeleton[row + 1, col] +
114
                                            64 * skeleton[row + 1, col - 1] +
115
                                            128 * skeleton[row, col - 1]]
116

117
                            if (neighbors == 0):
118
                                continue
119
                            elif ((neighbors == 3) or
120
                                  (neighbors == 1 and first_pass) or
121
                                  (neighbors == 2 and not first_pass)):
122
                                # Remove the pixel
123
                                cleaned_skeleton[row, col] = 0
124
                                pixel_removed = True
125

126
                # once a step has been processed, the original skeleton
127
                # is overwritten with the cleaned version
128
                skeleton[:, :] = cleaned_skeleton[:, :]
129

130
    return skeleton.base[1:-1, 1:-1].astype(bool)
131

132

133
def _skeletonize_loop(cnp.uint8_t[:, ::1] result,
134
                      Py_ssize_t[::1] i, Py_ssize_t[::1] j,
135
                      cnp.int32_t[::1] order, cnp.uint8_t[::1] table):
136
    """
137
    Inner loop of skeletonize function
138

139
    Parameters
140
    ----------
141

142
    result : ndarray of uint8
143
        On input, the image to be skeletonized, on output the skeletonized
144
        image.
145

146
    i, j : ndarrays
147
        The coordinates of each foreground pixel in the image
148

149
    order : ndarray
150
        The index of each pixel, in the order of processing (order[0] is
151
        the first pixel to process, etc.)
152

153
    table : ndarray
154
        The 512-element lookup table of values after transformation
155
        (whether to keep or not each configuration in a binary 3x3 array)
156

157
    Notes
158
    -----
159

160
    The loop determines whether each pixel in the image can be removed without
161
    changing the Euler number of the image. The pixels are ordered by
162
    increasing distance from the background which means a point nearer to
163
    the quench-line of the brushfire will be evaluated later than a
164
    point closer to the edge.
165

166
    Note that the neighborhood of a pixel may evolve before the loop
167
    arrives at this pixel. This is why it is possible to compute the
168
    skeleton in only one pass, thanks to an adapted ordering of the
169
    pixels.
170
    """
171
    cdef:
172
        Py_ssize_t accumulator
173
        Py_ssize_t index, order_index
174
        Py_ssize_t ii, jj
175
        Py_ssize_t rows = result.shape[0]
176
        Py_ssize_t cols = result.shape[1]
177

178
    with nogil:
179
        for index in range(order.shape[0]):
180
            accumulator = 16
181
            order_index = order[index]
182
            ii = i[order_index]
183
            jj = j[order_index]
184
            # Compute the configuration around the pixel
185
            if ii > 0:
186
                if jj > 0 and result[ii - 1, jj - 1]:
187
                    accumulator += 1
188
                if result[ii - 1, jj]:
189
                    accumulator += 2
190
                if jj < cols - 1 and result[ii - 1, jj + 1]:
191
                        accumulator += 4
192
            if jj > 0 and result[ii, jj - 1]:
193
                accumulator += 8
194
            if jj < cols - 1 and result[ii, jj + 1]:
195
                accumulator += 32
196
            if ii < rows - 1:
197
                if jj > 0 and result[ii + 1, jj - 1]:
198
                    accumulator += 64
199
                if result[ii + 1, jj]:
200
                    accumulator += 128
201
                if jj < cols - 1 and result[ii + 1, jj + 1]:
202
                    accumulator += 256
203
            # Assign the value of table corresponding to the configuration
204
            result[ii, jj] = table[accumulator]
205

206

207

208
def _table_lookup_index(cnp.uint8_t[:, ::1] image):
209
    """
210
    Return an index into a table per pixel of a binary image
211

212
    Take the sum of true neighborhood pixel values where the neighborhood
213
    looks like this::
214

215
       1   2   4
216
       8  16  32
217
      64 128 256
218

219
    This code could be replaced by a convolution with the kernel::
220

221
      256 128 64
222
       32  16  8
223
        4   2  1
224

225
    but this runs about twice as fast because of inlining and the
226
    hardwired kernel.
227
    """
228
    cdef:
229
        Py_ssize_t[:, ::1] indexer
230
        Py_ssize_t *p_indexer
231
        cnp.uint8_t *p_image
232
        Py_ssize_t i_stride
233
        Py_ssize_t i_shape
234
        Py_ssize_t j_shape
235
        Py_ssize_t i
236
        Py_ssize_t j
237
        Py_ssize_t offset
238

239
    i_shape   = image.shape[0]
240
    j_shape   = image.shape[1]
241
    indexer = np.zeros((i_shape, j_shape), dtype=np.intp)
242
    p_indexer = &indexer[0, 0]
243
    p_image   = &image[0, 0]
244
    i_stride  = image.strides[0]
245
    assert i_shape >= 3 and j_shape >= 3, \
246
        "Please use the slow method for arrays < 3x3"
247
    with nogil:
248
        for i in range(1, i_shape-1):
249
            offset = i_stride* i + 1
250
            for j in range(1, j_shape - 1):
251
                if p_image[offset]:
252
                    p_indexer[offset + i_stride + 1] += 1
253
                    p_indexer[offset + i_stride] += 2
254
                    p_indexer[offset + i_stride - 1] += 4
255
                    p_indexer[offset + 1] += 8
256
                    p_indexer[offset] += 16
257
                    p_indexer[offset - 1] += 32
258
                    p_indexer[offset - i_stride + 1] += 64
259
                    p_indexer[offset - i_stride] += 128
260
                    p_indexer[offset - i_stride - 1] += 256
261
                offset += 1
262
        #
263
        # Do the corner cases (literally)
264
        #
265
        if image[0, 0]:
266
            indexer[0, 0] += 16
267
            indexer[0, 1] += 8
268
            indexer[1, 0] += 2
269
            indexer[1, 1] += 1
270

271
        if image[0, j_shape - 1]:
272
            indexer[0, j_shape - 2] += 32
273
            indexer[0, j_shape - 1] += 16
274
            indexer[1, j_shape - 2] += 4
275
            indexer[1, j_shape - 1] += 2
276

277
        if image[i_shape - 1, 0]:
278
            indexer[i_shape - 2, 0] += 128
279
            indexer[i_shape - 2, 1] += 64
280
            indexer[i_shape - 1, 0] += 16
281
            indexer[i_shape - 1, 1] += 8
282

283
        if image[i_shape - 1, j_shape - 1]:
284
            indexer[i_shape - 2, j_shape - 2] += 256
285
            indexer[i_shape - 2, j_shape - 1] += 128
286
            indexer[i_shape - 1, j_shape - 2] += 32
287
            indexer[i_shape - 1, j_shape - 1] += 16
288
        #
289
        # Do the edges
290
        #
291
        for j in range(1, j_shape - 1):
292
            if image[0, j]:
293
                indexer[0, j - 1] += 32
294
                indexer[0, j] += 16
295
                indexer[0, j + 1] += 8
296
                indexer[1, j - 1] += 4
297
                indexer[1, j] += 2
298
                indexer[1, j + 1] += 1
299
            if image[i_shape - 1, j]:
300
                indexer[i_shape - 2, j - 1] += 256
301
                indexer[i_shape - 2, j] += 128
302
                indexer[i_shape - 2, j + 1] += 64
303
                indexer[i_shape - 1, j - 1] += 32
304
                indexer[i_shape - 1, j] += 16
305
                indexer[i_shape - 1, j + 1] += 8
306

307
        for i in range(1, i_shape - 1):
308
            if image[i, 0]:
309
                indexer[i - 1, 0] += 128
310
                indexer[i, 0] += 16
311
                indexer[i + 1, 0] += 2
312
                indexer[i - 1, 1] += 64
313
                indexer[i, 1] += 8
314
                indexer[i + 1, 1] += 1
315
            if image[i, j_shape - 1]:
316
                indexer[i - 1, j_shape - 2] += 256
317
                indexer[i, j_shape - 2] += 32
318
                indexer[i + 1, j_shape - 2] += 4
319
                indexer[i - 1, j_shape - 1] += 128
320
                indexer[i, j_shape - 1] += 16
321
                indexer[i + 1, j_shape - 1] += 2
322
    return np.asarray(indexer)
323