scikit-image

1
#cython: cdivision=True
2
#cython: boundscheck=False
3
#cython: nonecheck=False
4
#cython: wraparound=False
5
import numpy as np
6
cimport numpy as cnp
7
from libc.math cimport sin, cos
8
from .._shared.interpolation cimport bilinear_interpolation, round
9
from .._shared.transform cimport integrate
10

11
cdef extern from "numpy/npy_math.h":
12
    cnp.float64_t NAN "NPY_NAN"
13

14
from .._shared.fused_numerics cimport np_anyint as any_int
15

16
cnp.import_array()
17

18
def _glcm_loop(any_int[:, ::1] image, cnp.float64_t[:] distances,
19
               cnp.float64_t[:] angles, Py_ssize_t levels,
20
               cnp.uint32_t[:, :, :, ::1] out):
21
    """Perform co-occurrence matrix accumulation.
22

23
    Parameters
24
    ----------
25
    image : ndarray
26
        Integer typed input image. Only positive valued images are supported.
27
        If type is other than uint8, the argument `levels` needs to be set.
28
    distances : ndarray
29
        List of pixel pair distance offsets.
30
    angles : ndarray
31
        List of pixel pair angles in radians.
32
    levels : int
33
        The input image should contain integers in [0, `levels`-1],
34
        where levels indicate the number of gray-levels counted
35
        (typically 256 for an 8-bit image).
36
    out : ndarray
37
        On input a 4D array of zeros, and on output it contains
38
        the results of the GLCM computation.
39

40
    """
41

42
    cdef:
43
        Py_ssize_t a_idx, d_idx, r, c, rows, cols, row, col, start_row,\
44
                   end_row, start_col, end_col, offset_row, offset_col
45
        any_int i, j
46
        cnp.float64_t angle, distance
47

48
    with nogil:
49
        rows = image.shape[0]
50
        cols = image.shape[1]
51

52
        for a_idx in range(angles.shape[0]):
53
            angle = angles[a_idx]
54
            for d_idx in range(distances.shape[0]):
55
                distance = distances[d_idx]
56
                offset_row = round(sin(angle) * distance)
57
                offset_col = round(cos(angle) * distance)
58
                start_row = max(0, -offset_row)
59
                end_row = min(rows, rows - offset_row)
60
                start_col = max(0, -offset_col)
61
                end_col = min(cols, cols - offset_col)
62
                for r in range(start_row, end_row):
63
                    for c in range(start_col, end_col):
64
                        i = image[r, c]
65
                        # compute the location of the offset pixel
66
                        row = r + offset_row
67
                        col = c + offset_col
68
                        j = image[row, col]
69
                        if 0 <= i < levels and 0 <= j < levels:
70
                            out[i, j, d_idx, a_idx] += 1
71

72

73
cdef inline int _bit_rotate_right(int value, int length) noexcept nogil:
74
    """Cyclic bit shift to the right.
75

76
    Parameters
77
    ----------
78
    value : int
79
        integer value to shift
80
    length : int
81
        number of bits of integer
82

83
    """
84
    return (value >> 1) | ((value & 1) << (length - 1))
85

86

87
def _local_binary_pattern(cnp.float64_t[:, ::1] image,
88
                          int P, cnp.float64_t R, char method=b'D'):
89
    """Gray scale and rotation invariant LBP (Local Binary Patterns).
90

91
    LBP is an invariant descriptor that can be used for texture classification.
92

93
    Parameters
94
    ----------
95
    image : (N, M) cnp.float64_t array
96
        Graylevel image.
97
    P : int
98
        Number of circularly symmetric neighbor set points (quantization of
99
        the angular space).
100
    R : float
101
        Radius of circle (spatial resolution of the operator).
102
    method : {'D', 'R', 'U', 'N', 'V'}
103
        Method to determine the pattern.
104

105
        * 'D': 'default'
106
        * 'R': 'ror'
107
        * 'U': 'uniform'
108
        * 'N': 'nri_uniform'
109
        * 'V': 'var'
110

111
    Returns
112
    -------
113
    output : (N, M) array
114
        LBP image.
115
    """
116

117
    # texture weights
118
    cdef int[::1] weights = 2 ** np.arange(P, dtype=np.int32)
119
    # local position of texture elements
120
    rr = - R * np.sin(2 * np.pi * np.arange(P, dtype=np.float64) / P)
121
    cc = R * np.cos(2 * np.pi * np.arange(P, dtype=np.float64) / P)
122
    cdef cnp.float64_t[::1] rp = np.round(rr, 5)
123
    cdef cnp.float64_t[::1] cp = np.round(cc, 5)
124

125
    # pre-allocate arrays for computation
126
    cdef cnp.float64_t[::1] texture = np.zeros(P, dtype=np.float64)
127
    cdef signed char[::1] signed_texture = np.zeros(P, dtype=np.int8)
128
    cdef int[::1] rotation_chain = np.zeros(P, dtype=np.int32)
129

130
    output_shape = (image.shape[0], image.shape[1])
131
    cdef cnp.float64_t[:, ::1] output = np.zeros(output_shape, dtype=np.float64)
132

133
    cdef Py_ssize_t rows = image.shape[0]
134
    cdef Py_ssize_t cols = image.shape[1]
135

136
    cdef cnp.float64_t lbp
137
    cdef Py_ssize_t r, c, changes, i
138
    cdef Py_ssize_t rot_index, n_ones
139
    cdef cnp.int8_t first_zero, first_one
140

141
    # To compute the variance features
142
    cdef cnp.float64_t sum_, var_, texture_i
143

144
    with nogil:
145
        for r in range(image.shape[0]):
146
            for c in range(image.shape[1]):
147
                for i in range(P):
148
                    bilinear_interpolation[cnp.float64_t, cnp.float64_t, cnp.float64_t](
149
                            &image[0, 0], rows, cols, r + rp[i], c + cp[i],
150
                            b'C', 0, &texture[i])
151
                # signed / thresholded texture
152
                for i in range(P):
153
                    if texture[i] - image[r, c] >= 0:
154
                        signed_texture[i] = 1
155
                    else:
156
                        signed_texture[i] = 0
157

158
                lbp = 0
159

160
                # if method == b'var':
161
                if method == b'V':
162
                    # Compute the variance without passing from numpy.
163
                    # Following the LBP paper, we're taking a biased estimate
164
                    # of the variance (ddof=0)
165
                    sum_ = 0.0
166
                    var_ = 0.0
167
                    for i in range(P):
168
                        texture_i = texture[i]
169
                        sum_ += texture_i
170
                        var_ += texture_i * texture_i
171
                    var_ = (var_ - (sum_ * sum_) / P) / P
172
                    if var_ != 0:
173
                        lbp = var_
174
                    else:
175
                        lbp = NAN
176
                # if method == b'uniform':
177
                elif method == b'U' or method == b'N':
178
                    # determine number of 0 - 1 changes
179
                    changes = 0
180
                    for i in range(P - 1):
181
                        changes += (signed_texture[i]
182
                                    - signed_texture[i + 1]) != 0
183
                    if method == b'N':
184
                        # Uniform local binary patterns are defined as patterns
185
                        # with at most 2 value changes (from 0 to 1 or from 1 to
186
                        # 0). Uniform patterns can be characterized by their
187
                        # number `n_ones` of 1.  The possible values for
188
                        # `n_ones` range from 0 to P.
189
                        #
190
                        # Here is an example for P = 4:
191
                        # n_ones=0: 0000
192
                        # n_ones=1: 0001, 1000, 0100, 0010
193
                        # n_ones=2: 0011, 1001, 1100, 0110
194
                        # n_ones=3: 0111, 1011, 1101, 1110
195
                        # n_ones=4: 1111
196
                        #
197
                        # For a pattern of size P there are 2 constant patterns
198
                        # corresponding to n_ones=0 and n_ones=P. For each other
199
                        # value of `n_ones` , i.e n_ones=[1..P-1], there are P
200
                        # possible patterns which are related to each other
201
                        # through circular permutations. The total number of
202
                        # uniform patterns is thus (2 + P * (P - 1)).
203

204
                        # Given any pattern (uniform or not) we must be able to
205
                        # associate a unique code:
206
                        #
207
                        # 1. Constant patterns patterns (with n_ones=0 and
208
                        # n_ones=P) and non uniform patterns are given fixed
209
                        # code values.
210
                        #
211
                        # 2. Other uniform patterns are indexed considering the
212
                        # value of n_ones, and an index called 'rot_index'
213
                        # representing the number of circular right shifts
214
                        # required to obtain the pattern starting from a
215
                        # reference position (corresponding to all zeros stacked
216
                        # on the right). This number of rotations (or circular
217
                        # right shifts) 'rot_index' is efficiently computed by
218
                        # considering the positions of the first 1 and the first
219
                        # 0 found in the pattern.
220

221
                        if changes <= 2:
222
                            # We have a uniform pattern
223
                            n_ones = 0  # determines the number of ones
224
                            first_one = -1  # position was the first one
225
                            first_zero = -1  # position of the first zero
226
                            for i in range(P):
227
                                if signed_texture[i]:
228
                                    n_ones += 1
229
                                    if first_one == -1:
230
                                        first_one = i
231
                                else:
232
                                    if first_zero == -1:
233
                                        first_zero = i
234
                            if n_ones == 0:
235
                                lbp = 0
236
                            elif n_ones == P:
237
                                lbp = P * (P - 1) + 1
238
                            else:
239
                                if first_one == 0:
240
                                    rot_index = n_ones - first_zero
241
                                else:
242
                                    rot_index = P - first_one
243
                                lbp = 1 + (n_ones - 1) * P + rot_index
244
                        else:  # changes > 2
245
                            lbp = P * (P - 1) + 2
246
                    else:  # method != 'N'
247
                        if changes <= 2:
248
                            for i in range(P):
249
                                lbp += signed_texture[i]
250
                        else:
251
                            lbp = P + 1
252
                else:
253
                    # method == b'default'
254
                    for i in range(P):
255
                        lbp += signed_texture[i] * weights[i]
256

257
                    # method == b'ror'
258
                    if method == b'R':
259
                        # shift LBP P times to the right and get minimum value
260
                        rotation_chain[0] = <int>lbp
261
                        for i in range(1, P):
262
                            rotation_chain[i] = \
263
                                _bit_rotate_right(rotation_chain[i - 1], P)
264
                        lbp = rotation_chain[0]
265
                        for i in range(1, P):
266
                            lbp = min(lbp, rotation_chain[i])
267

268
                output[r, c] = lbp
269

270
    return np.asarray(output)
271

272

273
# Constant values that are used by `_multiblock_lbp` function.
274
# Values represent offsets of neighbor rectangles relative to central one.
275
# It has order starting from top left and going clockwise.
276
cdef:
277
    Py_ssize_t[::1] mlbp_r_offsets = np.asarray([-1, -1, -1, 0, 1, 1, 1, 0],
278
                                                dtype=np.intp)
279
    Py_ssize_t[::1] mlbp_c_offsets = np.asarray([-1, 0, 1, 1, 1, 0, -1, -1],
280
                                                dtype=np.intp)
281

282

283
cpdef int _multiblock_lbp(np_floats[:, ::1] int_image,
284
                          Py_ssize_t r,
285
                          Py_ssize_t c,
286
                          Py_ssize_t width,
287
                          Py_ssize_t height) noexcept nogil:
288
    """Multi-block local binary pattern (MB-LBP) [1]_.
289

290
    Parameters
291
    ----------
292
    int_image : (N, M) float array
293
        Integral image.
294
    r : int
295
        Row-coordinate of top left corner of a rectangle containing feature.
296
    c : int
297
        Column-coordinate of top left corner of a rectangle containing feature.
298
    width : int
299
        Width of one of 9 equal rectangles that will be used to compute
300
        a feature.
301
    height : int
302
        Height of one of 9 equal rectangles that will be used to compute
303
        a feature.
304

305
    Returns
306
    -------
307
    output : int
308
        8-bit MB-LBP feature descriptor.
309

310
    References
311
    ----------
312
    .. [1] L. Zhang, R. Chu, S. Xiang, S. Liao, S.Z. Li. "Face Detection Based
313
           on Multi-Block LBP Representation", In Proceedings: Advances in
314
           Biometrics, International Conference, ICB 2007, Seoul, Korea.
315
           http://www.cbsr.ia.ac.cn/users/scliao/papers/Zhang-ICB07-MBLBP.pdf
316
           :DOI:`10.1007/978-3-540-74549-5_2`
317
    """
318

319
    cdef:
320
        # Top-left coordinates of central rectangle.
321
        Py_ssize_t central_rect_r = r + height
322
        Py_ssize_t central_rect_c = c + width
323

324
        Py_ssize_t r_shift = height - 1
325
        Py_ssize_t c_shift = width - 1
326

327
        Py_ssize_t current_rect_r, current_rect_c
328
        Py_ssize_t element_num
329
        np_floats current_rect_val
330
        int has_greater_value
331
        int lbp_code = 0
332

333
    # Sum of intensity values of central rectangle.
334
    cdef np_floats central_rect_val = integrate(int_image, central_rect_r,
335
                                                central_rect_c,
336
                                                central_rect_r + r_shift,
337
                                                central_rect_c + c_shift)
338

339
    for element_num in range(8):
340

341
        current_rect_r = central_rect_r + mlbp_r_offsets[element_num]*height
342
        current_rect_c = central_rect_c + mlbp_c_offsets[element_num]*width
343

344

345
        current_rect_val = integrate(int_image, current_rect_r, current_rect_c,
346
                                     current_rect_r + r_shift,
347
                                     current_rect_c + c_shift)
348

349

350
        has_greater_value = current_rect_val >= central_rect_val
351

352
        # If current rectangle's intensity value is bigger
353
        # make corresponding bit to 1.
354
        lbp_code |= has_greater_value << (7 - element_num)
355

356
    return lbp_code
357