-
-
Notifications
You must be signed in to change notification settings - Fork 2.2k
/
_texture.pyx
356 lines (304 loc) · 13.9 KB
/
_texture.pyx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as cnp
from libc.math cimport sin, cos
from .._shared.interpolation cimport bilinear_interpolation, round
from .._shared.transform cimport integrate
cdef extern from "numpy/npy_math.h":
cnp.float64_t NAN "NPY_NAN"
from .._shared.fused_numerics cimport np_anyint as any_int
cnp.import_array()
def _glcm_loop(any_int[:, ::1] image, cnp.float64_t[:] distances,
cnp.float64_t[:] angles, Py_ssize_t levels,
cnp.uint32_t[:, :, :, ::1] out):
"""Perform co-occurrence matrix accumulation.
Parameters
----------
image : ndarray
Integer typed input image. Only positive valued images are supported.
If type is other than uint8, the argument `levels` needs to be set.
distances : ndarray
List of pixel pair distance offsets.
angles : ndarray
List of pixel pair angles in radians.
levels : int
The input image should contain integers in [0, `levels`-1],
where levels indicate the number of gray-levels counted
(typically 256 for an 8-bit image).
out : ndarray
On input a 4D array of zeros, and on output it contains
the results of the GLCM computation.
"""
cdef:
Py_ssize_t a_idx, d_idx, r, c, rows, cols, row, col, start_row,\
end_row, start_col, end_col, offset_row, offset_col
any_int i, j
cnp.float64_t angle, distance
with nogil:
rows = image.shape[0]
cols = image.shape[1]
for a_idx in range(angles.shape[0]):
angle = angles[a_idx]
for d_idx in range(distances.shape[0]):
distance = distances[d_idx]
offset_row = round(sin(angle) * distance)
offset_col = round(cos(angle) * distance)
start_row = max(0, -offset_row)
end_row = min(rows, rows - offset_row)
start_col = max(0, -offset_col)
end_col = min(cols, cols - offset_col)
for r in range(start_row, end_row):
for c in range(start_col, end_col):
i = image[r, c]
# compute the location of the offset pixel
row = r + offset_row
col = c + offset_col
j = image[row, col]
if 0 <= i < levels and 0 <= j < levels:
out[i, j, d_idx, a_idx] += 1
cdef inline int _bit_rotate_right(int value, int length) noexcept nogil:
"""Cyclic bit shift to the right.
Parameters
----------
value : int
integer value to shift
length : int
number of bits of integer
"""
return (value >> 1) | ((value & 1) << (length - 1))
def _local_binary_pattern(cnp.float64_t[:, ::1] image,
int P, cnp.float64_t R, char method=b'D'):
"""Gray scale and rotation invariant LBP (Local Binary Patterns).
LBP is an invariant descriptor that can be used for texture classification.
Parameters
----------
image : (N, M) cnp.float64_t array
Graylevel image.
P : int
Number of circularly symmetric neighbor set points (quantization of
the angular space).
R : float
Radius of circle (spatial resolution of the operator).
method : {'D', 'R', 'U', 'N', 'V'}
Method to determine the pattern.
* 'D': 'default'
* 'R': 'ror'
* 'U': 'uniform'
* 'N': 'nri_uniform'
* 'V': 'var'
Returns
-------
output : (N, M) array
LBP image.
"""
# texture weights
cdef int[::1] weights = 2 ** np.arange(P, dtype=np.int32)
# local position of texture elements
rr = - R * np.sin(2 * np.pi * np.arange(P, dtype=np.float64) / P)
cc = R * np.cos(2 * np.pi * np.arange(P, dtype=np.float64) / P)
cdef cnp.float64_t[::1] rp = np.round(rr, 5)
cdef cnp.float64_t[::1] cp = np.round(cc, 5)
# pre-allocate arrays for computation
cdef cnp.float64_t[::1] texture = np.zeros(P, dtype=np.float64)
cdef signed char[::1] signed_texture = np.zeros(P, dtype=np.int8)
cdef int[::1] rotation_chain = np.zeros(P, dtype=np.int32)
output_shape = (image.shape[0], image.shape[1])
cdef cnp.float64_t[:, ::1] output = np.zeros(output_shape, dtype=np.float64)
cdef Py_ssize_t rows = image.shape[0]
cdef Py_ssize_t cols = image.shape[1]
cdef cnp.float64_t lbp
cdef Py_ssize_t r, c, changes, i
cdef Py_ssize_t rot_index, n_ones
cdef cnp.int8_t first_zero, first_one
# To compute the variance features
cdef cnp.float64_t sum_, var_, texture_i
with nogil:
for r in range(image.shape[0]):
for c in range(image.shape[1]):
for i in range(P):
bilinear_interpolation[cnp.float64_t, cnp.float64_t, cnp.float64_t](
&image[0, 0], rows, cols, r + rp[i], c + cp[i],
b'C', 0, &texture[i])
# signed / thresholded texture
for i in range(P):
if texture[i] - image[r, c] >= 0:
signed_texture[i] = 1
else:
signed_texture[i] = 0
lbp = 0
# if method == b'var':
if method == b'V':
# Compute the variance without passing from numpy.
# Following the LBP paper, we're taking a biased estimate
# of the variance (ddof=0)
sum_ = 0.0
var_ = 0.0
for i in range(P):
texture_i = texture[i]
sum_ += texture_i
var_ += texture_i * texture_i
var_ = (var_ - (sum_ * sum_) / P) / P
if var_ != 0:
lbp = var_
else:
lbp = NAN
# if method == b'uniform':
elif method == b'U' or method == b'N':
# determine number of 0 - 1 changes
changes = 0
for i in range(P - 1):
changes += (signed_texture[i]
- signed_texture[i + 1]) != 0
if method == b'N':
# Uniform local binary patterns are defined as patterns
# with at most 2 value changes (from 0 to 1 or from 1 to
# 0). Uniform patterns can be characterized by their
# number `n_ones` of 1. The possible values for
# `n_ones` range from 0 to P.
#
# Here is an example for P = 4:
# n_ones=0: 0000
# n_ones=1: 0001, 1000, 0100, 0010
# n_ones=2: 0011, 1001, 1100, 0110
# n_ones=3: 0111, 1011, 1101, 1110
# n_ones=4: 1111
#
# For a pattern of size P there are 2 constant patterns
# corresponding to n_ones=0 and n_ones=P. For each other
# value of `n_ones` , i.e n_ones=[1..P-1], there are P
# possible patterns which are related to each other
# through circular permutations. The total number of
# uniform patterns is thus (2 + P * (P - 1)).
# Given any pattern (uniform or not) we must be able to
# associate a unique code:
#
# 1. Constant patterns patterns (with n_ones=0 and
# n_ones=P) and non uniform patterns are given fixed
# code values.
#
# 2. Other uniform patterns are indexed considering the
# value of n_ones, and an index called 'rot_index'
# representing the number of circular right shifts
# required to obtain the pattern starting from a
# reference position (corresponding to all zeros stacked
# on the right). This number of rotations (or circular
# right shifts) 'rot_index' is efficiently computed by
# considering the positions of the first 1 and the first
# 0 found in the pattern.
if changes <= 2:
# We have a uniform pattern
n_ones = 0 # determines the number of ones
first_one = -1 # position was the first one
first_zero = -1 # position of the first zero
for i in range(P):
if signed_texture[i]:
n_ones += 1
if first_one == -1:
first_one = i
else:
if first_zero == -1:
first_zero = i
if n_ones == 0:
lbp = 0
elif n_ones == P:
lbp = P * (P - 1) + 1
else:
if first_one == 0:
rot_index = n_ones - first_zero
else:
rot_index = P - first_one
lbp = 1 + (n_ones - 1) * P + rot_index
else: # changes > 2
lbp = P * (P - 1) + 2
else: # method != 'N'
if changes <= 2:
for i in range(P):
lbp += signed_texture[i]
else:
lbp = P + 1
else:
# method == b'default'
for i in range(P):
lbp += signed_texture[i] * weights[i]
# method == b'ror'
if method == b'R':
# shift LBP P times to the right and get minimum value
rotation_chain[0] = <int>lbp
for i in range(1, P):
rotation_chain[i] = \
_bit_rotate_right(rotation_chain[i - 1], P)
lbp = rotation_chain[0]
for i in range(1, P):
lbp = min(lbp, rotation_chain[i])
output[r, c] = lbp
return np.asarray(output)
# Constant values that are used by `_multiblock_lbp` function.
# Values represent offsets of neighbor rectangles relative to central one.
# It has order starting from top left and going clockwise.
cdef:
Py_ssize_t[::1] mlbp_r_offsets = np.asarray([-1, -1, -1, 0, 1, 1, 1, 0],
dtype=np.intp)
Py_ssize_t[::1] mlbp_c_offsets = np.asarray([-1, 0, 1, 1, 1, 0, -1, -1],
dtype=np.intp)
cpdef int _multiblock_lbp(np_floats[:, ::1] int_image,
Py_ssize_t r,
Py_ssize_t c,
Py_ssize_t width,
Py_ssize_t height) noexcept nogil:
"""Multi-block local binary pattern (MB-LBP) [1]_.
Parameters
----------
int_image : (N, M) float array
Integral image.
r : int
Row-coordinate of top left corner of a rectangle containing feature.
c : int
Column-coordinate of top left corner of a rectangle containing feature.
width : int
Width of one of 9 equal rectangles that will be used to compute
a feature.
height : int
Height of one of 9 equal rectangles that will be used to compute
a feature.
Returns
-------
output : int
8-bit MB-LBP feature descriptor.
References
----------
.. [1] L. Zhang, R. Chu, S. Xiang, S. Liao, S.Z. Li. "Face Detection Based
on Multi-Block LBP Representation", In Proceedings: Advances in
Biometrics, International Conference, ICB 2007, Seoul, Korea.
http://www.cbsr.ia.ac.cn/users/scliao/papers/Zhang-ICB07-MBLBP.pdf
:DOI:`10.1007/978-3-540-74549-5_2`
"""
cdef:
# Top-left coordinates of central rectangle.
Py_ssize_t central_rect_r = r + height
Py_ssize_t central_rect_c = c + width
Py_ssize_t r_shift = height - 1
Py_ssize_t c_shift = width - 1
Py_ssize_t current_rect_r, current_rect_c
Py_ssize_t element_num
np_floats current_rect_val
int has_greater_value
int lbp_code = 0
# Sum of intensity values of central rectangle.
cdef np_floats central_rect_val = integrate(int_image, central_rect_r,
central_rect_c,
central_rect_r + r_shift,
central_rect_c + c_shift)
for element_num in range(8):
current_rect_r = central_rect_r + mlbp_r_offsets[element_num]*height
current_rect_c = central_rect_c + mlbp_c_offsets[element_num]*width
current_rect_val = integrate(int_image, current_rect_r, current_rect_c,
current_rect_r + r_shift,
current_rect_c + c_shift)
has_greater_value = current_rect_val >= central_rect_val
# If current rectangle's intensity value is bigger
# make corresponding bit to 1.
lbp_code |= has_greater_value << (7 - element_num)
return lbp_code