texture.py 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562
  1. """
  2. Methods to characterize image textures.
  3. """
  4. import warnings
  5. import numpy as np
  6. from .._shared.utils import check_nD
  7. from ..color import gray2rgb
  8. from ..util import img_as_float
  9. from ._texture import _glcm_loop, _local_binary_pattern, _multiblock_lbp
  10. def graycomatrix(image, distances, angles, levels=None, symmetric=False, normed=False):
  11. """Calculate the gray-level co-occurrence matrix.
  12. A gray level co-occurrence matrix is a histogram of co-occurring
  13. grayscale values at a given offset over an image.
  14. .. versionchanged:: 0.19
  15. `greymatrix` was renamed to `graymatrix` in 0.19.
  16. Parameters
  17. ----------
  18. image : array_like
  19. Integer typed input image. Only positive valued images are supported.
  20. If type is other than uint8, the argument `levels` needs to be set.
  21. distances : array_like
  22. List of pixel pair distance offsets.
  23. angles : array_like
  24. List of pixel pair angles in radians.
  25. levels : int, optional
  26. The input image should contain integers in [0, `levels`-1],
  27. where levels indicate the number of gray-levels counted
  28. (typically 256 for an 8-bit image). This argument is required for
  29. 16-bit images or higher and is typically the maximum of the image.
  30. As the output matrix is at least `levels` x `levels`, it might
  31. be preferable to use binning of the input image rather than
  32. large values for `levels`.
  33. symmetric : bool, optional
  34. If True, the output matrix `P[:, :, d, theta]` is symmetric. This
  35. is accomplished by ignoring the order of value pairs, so both
  36. (i, j) and (j, i) are accumulated when (i, j) is encountered
  37. for a given offset. The default is False.
  38. normed : bool, optional
  39. If True, normalize each matrix `P[:, :, d, theta]` by dividing
  40. by the total number of accumulated co-occurrences for the given
  41. offset. The elements of the resulting matrix sum to 1. The
  42. default is False.
  43. Returns
  44. -------
  45. P : 4-D ndarray
  46. The gray-level co-occurrence histogram. The value
  47. `P[i,j,d,theta]` is the number of times that gray-level `j`
  48. occurs at a distance `d` and at an angle `theta` from
  49. gray-level `i`. If `normed` is `False`, the output is of
  50. type uint32, otherwise it is float64. The dimensions are:
  51. levels x levels x number of distances x number of angles.
  52. References
  53. ----------
  54. .. [1] M. Hall-Beyer, 2007. GLCM Texture: A Tutorial
  55. https://prism.ucalgary.ca/handle/1880/51900
  56. DOI:`10.11575/PRISM/33280`
  57. .. [2] R.M. Haralick, K. Shanmugam, and I. Dinstein, "Textural features for
  58. image classification", IEEE Transactions on Systems, Man, and
  59. Cybernetics, vol. SMC-3, no. 6, pp. 610-621, Nov. 1973.
  60. :DOI:`10.1109/TSMC.1973.4309314`
  61. .. [3] M. Nadler and E.P. Smith, Pattern Recognition Engineering,
  62. Wiley-Interscience, 1993.
  63. .. [4] Wikipedia, https://en.wikipedia.org/wiki/Co-occurrence_matrix
  64. Examples
  65. --------
  66. Compute 4 GLCMs using 1-pixel distance and 4 different angles. For example,
  67. an angle of 0 radians refers to the neighboring pixel to the right;
  68. pi/4 radians to the top-right diagonal neighbor; pi/2 radians to the pixel
  69. above, and so forth.
  70. >>> image = np.array([[0, 0, 1, 1],
  71. ... [0, 0, 1, 1],
  72. ... [0, 2, 2, 2],
  73. ... [2, 2, 3, 3]], dtype=np.uint8)
  74. >>> result = graycomatrix(image, [1], [0, np.pi/4, np.pi/2, 3*np.pi/4],
  75. ... levels=4)
  76. >>> result[:, :, 0, 0]
  77. array([[2, 2, 1, 0],
  78. [0, 2, 0, 0],
  79. [0, 0, 3, 1],
  80. [0, 0, 0, 1]], dtype=uint32)
  81. >>> result[:, :, 0, 1]
  82. array([[1, 1, 3, 0],
  83. [0, 1, 1, 0],
  84. [0, 0, 0, 2],
  85. [0, 0, 0, 0]], dtype=uint32)
  86. >>> result[:, :, 0, 2]
  87. array([[3, 0, 2, 0],
  88. [0, 2, 2, 0],
  89. [0, 0, 1, 2],
  90. [0, 0, 0, 0]], dtype=uint32)
  91. >>> result[:, :, 0, 3]
  92. array([[2, 0, 0, 0],
  93. [1, 1, 2, 0],
  94. [0, 0, 2, 1],
  95. [0, 0, 0, 0]], dtype=uint32)
  96. """
  97. check_nD(image, 2)
  98. check_nD(distances, 1, 'distances')
  99. check_nD(angles, 1, 'angles')
  100. image = np.ascontiguousarray(image)
  101. image_max = image.max()
  102. if np.issubdtype(image.dtype, np.floating):
  103. raise ValueError(
  104. "Float images are not supported by graycomatrix. "
  105. "Convert the image to an unsigned integer type."
  106. )
  107. # for image type > 8bit, levels must be set.
  108. if image.dtype not in (np.uint8, np.int8) and levels is None:
  109. raise ValueError(
  110. "The levels argument is required for data types "
  111. "other than uint8. The resulting matrix will be at "
  112. "least levels ** 2 in size."
  113. )
  114. if np.issubdtype(image.dtype, np.signedinteger) and np.any(image < 0):
  115. raise ValueError("Negative-valued images are not supported.")
  116. if levels is None:
  117. levels = 256
  118. if image_max >= levels:
  119. raise ValueError(
  120. "The maximum grayscale value in the image should be "
  121. "smaller than the number of levels."
  122. )
  123. distances = np.ascontiguousarray(distances, dtype=np.float64)
  124. angles = np.ascontiguousarray(angles, dtype=np.float64)
  125. P = np.zeros(
  126. (levels, levels, len(distances), len(angles)), dtype=np.uint32, order='C'
  127. )
  128. # count co-occurences
  129. _glcm_loop(image, distances, angles, levels, P)
  130. # make each GLMC symmetric
  131. if symmetric:
  132. Pt = np.transpose(P, (1, 0, 2, 3))
  133. P = P + Pt
  134. # normalize each GLCM
  135. if normed:
  136. P = P.astype(np.float64)
  137. glcm_sums = np.sum(P, axis=(0, 1), keepdims=True)
  138. glcm_sums[glcm_sums == 0] = 1
  139. P /= glcm_sums
  140. return P
  141. def graycoprops(P, prop='contrast'):
  142. """Calculate texture properties of a GLCM.
  143. Compute a feature of a gray level co-occurrence matrix to serve as
  144. a compact summary of the matrix. The properties are computed as
  145. follows:
  146. - 'contrast': :math:`\\sum_{i,j=0}^{levels-1} P_{i,j}(i-j)^2`
  147. - 'dissimilarity': :math:`\\sum_{i,j=0}^{levels-1}P_{i,j}|i-j|`
  148. - 'homogeneity': :math:`\\sum_{i,j=0}^{levels-1}\\frac{P_{i,j}}{1+(i-j)^2}`
  149. - 'ASM': :math:`\\sum_{i,j=0}^{levels-1} P_{i,j}^2`
  150. - 'energy': :math:`\\sqrt{ASM}`
  151. - 'correlation':
  152. .. math:: \\sum_{i,j=0}^{levels-1} P_{i,j}\\left[\\frac{(i-\\mu_i) \\
  153. (j-\\mu_j)}{\\sqrt{(\\sigma_i^2)(\\sigma_j^2)}}\\right]
  154. - 'mean': :math:`\\sum_{i=0}^{levels-1} i*P_{i}`
  155. - 'variance': :math:`\\sum_{i=0}^{levels-1} P_{i}*(i-mean)^2`
  156. - 'std': :math:`\\sqrt{variance}`
  157. - 'entropy': :math:`\\sum_{i,j=0}^{levels-1} -P_{i,j}*log(P_{i,j})`
  158. Each GLCM is normalized to have a sum of 1 before the computation of
  159. texture properties.
  160. .. versionchanged:: 0.19
  161. `greycoprops` was renamed to `graycoprops` in 0.19.
  162. Parameters
  163. ----------
  164. P : ndarray
  165. Input array. `P` is the gray-level co-occurrence histogram
  166. for which to compute the specified property. The value
  167. `P[i,j,d,theta]` is the number of times that gray-level j
  168. occurs at a distance d and at an angle theta from
  169. gray-level i.
  170. prop : {'contrast', 'dissimilarity', 'homogeneity', 'energy', \
  171. 'correlation', 'ASM', 'mean', 'variance', 'std', 'entropy'}, optional
  172. The property of the GLCM to compute. The default is 'contrast'.
  173. Returns
  174. -------
  175. results : 2-D ndarray
  176. 2-dimensional array. `results[d, a]` is the property 'prop' for
  177. the d'th distance and the a'th angle.
  178. References
  179. ----------
  180. .. [1] M. Hall-Beyer, 2007. GLCM Texture: A Tutorial v. 1.0 through 3.0.
  181. The GLCM Tutorial Home Page,
  182. https://prism.ucalgary.ca/handle/1880/51900
  183. DOI:`10.11575/PRISM/33280`
  184. Examples
  185. --------
  186. Compute the contrast for GLCMs with distances [1, 2] and angles
  187. [0 degrees, 90 degrees]
  188. >>> image = np.array([[0, 0, 1, 1],
  189. ... [0, 0, 1, 1],
  190. ... [0, 2, 2, 2],
  191. ... [2, 2, 3, 3]], dtype=np.uint8)
  192. >>> g = graycomatrix(image, [1, 2], [0, np.pi/2], levels=4,
  193. ... normed=True, symmetric=True)
  194. >>> contrast = graycoprops(g, 'contrast')
  195. >>> contrast
  196. array([[0.58333333, 1. ],
  197. [1.25 , 2.75 ]])
  198. """
  199. def glcm_mean():
  200. I = np.arange(num_level).reshape((num_level, 1, 1, 1))
  201. mean = np.sum(I * P, axis=(0, 1))
  202. return I, mean
  203. check_nD(P, 4, 'P')
  204. (num_level, num_level2, num_dist, num_angle) = P.shape
  205. if num_level != num_level2:
  206. raise ValueError('num_level and num_level2 must be equal.')
  207. if num_dist <= 0:
  208. raise ValueError('num_dist must be positive.')
  209. if num_angle <= 0:
  210. raise ValueError('num_angle must be positive.')
  211. # normalize each GLCM
  212. P = P.astype(np.float64)
  213. glcm_sums = np.sum(P, axis=(0, 1), keepdims=True)
  214. glcm_sums[glcm_sums == 0] = 1
  215. P /= glcm_sums
  216. # create weights for specified property
  217. I, J = np.ogrid[0:num_level, 0:num_level]
  218. if prop == 'contrast':
  219. weights = (I - J) ** 2
  220. elif prop == 'dissimilarity':
  221. weights = np.abs(I - J)
  222. elif prop == 'homogeneity':
  223. weights = 1.0 / (1.0 + (I - J) ** 2)
  224. elif prop in ['ASM', 'energy', 'correlation', 'entropy', 'variance', 'mean', 'std']:
  225. pass
  226. else:
  227. raise ValueError(f'{prop} is an invalid property')
  228. # compute property for each GLCM
  229. if prop == 'energy':
  230. asm = np.sum(P**2, axis=(0, 1))
  231. results = np.sqrt(asm)
  232. elif prop == 'ASM':
  233. results = np.sum(P**2, axis=(0, 1))
  234. elif prop == 'mean':
  235. _, results = glcm_mean()
  236. elif prop == 'variance':
  237. I, mean = glcm_mean()
  238. results = np.sum(P * ((I - mean) ** 2), axis=(0, 1))
  239. elif prop == 'std':
  240. I, mean = glcm_mean()
  241. var = np.sum(P * ((I - mean) ** 2), axis=(0, 1))
  242. results = np.sqrt(var)
  243. elif prop == 'entropy':
  244. ln = -np.log(P, where=(P != 0), out=np.zeros_like(P))
  245. results = np.sum(P * ln, axis=(0, 1))
  246. elif prop == 'correlation':
  247. results = np.zeros((num_dist, num_angle), dtype=np.float64)
  248. I = np.array(range(num_level)).reshape((num_level, 1, 1, 1))
  249. J = np.array(range(num_level)).reshape((1, num_level, 1, 1))
  250. diff_i = I - np.sum(I * P, axis=(0, 1))
  251. diff_j = J - np.sum(J * P, axis=(0, 1))
  252. std_i = np.sqrt(np.sum(P * (diff_i) ** 2, axis=(0, 1)))
  253. std_j = np.sqrt(np.sum(P * (diff_j) ** 2, axis=(0, 1)))
  254. cov = np.sum(P * (diff_i * diff_j), axis=(0, 1))
  255. # handle the special case of standard deviations near zero
  256. mask_0 = std_i < 1e-15
  257. mask_0[std_j < 1e-15] = True
  258. results[mask_0] = 1
  259. # handle the standard case
  260. mask_1 = ~mask_0
  261. results[mask_1] = cov[mask_1] / (std_i[mask_1] * std_j[mask_1])
  262. elif prop in ['contrast', 'dissimilarity', 'homogeneity']:
  263. weights = weights.reshape((num_level, num_level, 1, 1))
  264. results = np.sum(P * weights, axis=(0, 1))
  265. return results
  266. def local_binary_pattern(image, P, R, method='default'):
  267. """Compute the local binary patterns (LBP) of an image.
  268. LBP is a visual descriptor often used in texture classification.
  269. Parameters
  270. ----------
  271. image : (M, N) array
  272. 2D grayscale image.
  273. P : int
  274. Number of circularly symmetric neighbor set points (quantization of
  275. the angular space).
  276. R : float
  277. Radius of circle (spatial resolution of the operator).
  278. method : str {'default', 'ror', 'uniform', 'nri_uniform', 'var'}, optional
  279. Method to determine the pattern:
  280. ``default``
  281. Original local binary pattern which is grayscale invariant but not
  282. rotation invariant.
  283. ``ror``
  284. Extension of default pattern which is grayscale invariant and
  285. rotation invariant.
  286. ``uniform``
  287. Uniform pattern which is grayscale invariant and rotation
  288. invariant, offering finer quantization of the angular space.
  289. For details, see [1]_.
  290. ``nri_uniform``
  291. Variant of uniform pattern which is grayscale invariant but not
  292. rotation invariant. For details, see [2]_ and [3]_.
  293. ``var``
  294. Variance of local image texture (related to contrast)
  295. which is rotation invariant but not grayscale invariant.
  296. Returns
  297. -------
  298. output : (M, N) array
  299. LBP image.
  300. References
  301. ----------
  302. .. [1] T. Ojala, M. Pietikainen, T. Maenpaa, "Multiresolution gray-scale
  303. and rotation invariant texture classification with local binary
  304. patterns", IEEE Transactions on Pattern Analysis and Machine
  305. Intelligence, vol. 24, no. 7, pp. 971-987, July 2002
  306. :DOI:`10.1109/TPAMI.2002.1017623`
  307. .. [2] T. Ahonen, A. Hadid and M. Pietikainen. "Face recognition with
  308. local binary patterns", in Proc. Eighth European Conf. Computer
  309. Vision, Prague, Czech Republic, May 11-14, 2004, pp. 469-481, 2004.
  310. http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.214.6851
  311. :DOI:`10.1007/978-3-540-24670-1_36`
  312. .. [3] T. Ahonen, A. Hadid and M. Pietikainen, "Face Description with
  313. Local Binary Patterns: Application to Face Recognition",
  314. IEEE Transactions on Pattern Analysis and Machine Intelligence,
  315. vol. 28, no. 12, pp. 2037-2041, Dec. 2006
  316. :DOI:`10.1109/TPAMI.2006.244`
  317. """
  318. check_nD(image, 2)
  319. methods = {
  320. 'default': ord('D'),
  321. 'ror': ord('R'),
  322. 'uniform': ord('U'),
  323. 'nri_uniform': ord('N'),
  324. 'var': ord('V'),
  325. }
  326. if np.issubdtype(image.dtype, np.floating):
  327. warnings.warn(
  328. "Applying `local_binary_pattern` to floating-point images may "
  329. "give unexpected results when small numerical differences between "
  330. "adjacent pixels are present. It is recommended to use this "
  331. "function with images of integer dtype."
  332. )
  333. image = np.ascontiguousarray(image, dtype=np.float64)
  334. output = _local_binary_pattern(image, P, R, methods[method.lower()])
  335. return output
  336. def multiblock_lbp(int_image, r, c, width, height):
  337. """Multi-block local binary pattern (MB-LBP).
  338. The features are calculated similarly to local binary patterns (LBPs),
  339. (See :py:meth:`local_binary_pattern`) except that summed blocks are
  340. used instead of individual pixel values.
  341. MB-LBP is an extension of LBP that can be computed on multiple scales
  342. in constant time using the integral image. Nine equally-sized rectangles
  343. are used to compute a feature. For each rectangle, the sum of the pixel
  344. intensities is computed. Comparisons of these sums to that of the central
  345. rectangle determine the feature, similarly to LBP.
  346. Parameters
  347. ----------
  348. int_image : (N, M) array
  349. Integral image.
  350. r : int
  351. Row-coordinate of top left corner of a rectangle containing feature.
  352. c : int
  353. Column-coordinate of top left corner of a rectangle containing feature.
  354. width : int
  355. Width of one of the 9 equal rectangles that will be used to compute
  356. a feature.
  357. height : int
  358. Height of one of the 9 equal rectangles that will be used to compute
  359. a feature.
  360. Returns
  361. -------
  362. output : int
  363. 8-bit MB-LBP feature descriptor.
  364. References
  365. ----------
  366. .. [1] L. Zhang, R. Chu, S. Xiang, S. Liao, S.Z. Li. "Face Detection Based
  367. on Multi-Block LBP Representation", In Proceedings: Advances in
  368. Biometrics, International Conference, ICB 2007, Seoul, Korea.
  369. http://www.cbsr.ia.ac.cn/users/scliao/papers/Zhang-ICB07-MBLBP.pdf
  370. :DOI:`10.1007/978-3-540-74549-5_2`
  371. """
  372. int_image = np.ascontiguousarray(int_image, dtype=np.float32)
  373. lbp_code = _multiblock_lbp(int_image, r, c, width, height)
  374. return lbp_code
  375. def draw_multiblock_lbp(
  376. image,
  377. r,
  378. c,
  379. width,
  380. height,
  381. lbp_code=0,
  382. color_greater_block=(1, 1, 1),
  383. color_less_block=(0, 0.69, 0.96),
  384. alpha=0.5,
  385. ):
  386. """Multi-block local binary pattern visualization.
  387. Blocks with higher sums are colored with alpha-blended white rectangles,
  388. whereas blocks with lower sums are colored alpha-blended cyan. Colors
  389. and the `alpha` parameter can be changed.
  390. Parameters
  391. ----------
  392. image : ndarray of float or uint
  393. Image on which to visualize the pattern.
  394. r : int
  395. Row-coordinate of top left corner of a rectangle containing feature.
  396. c : int
  397. Column-coordinate of top left corner of a rectangle containing feature.
  398. width : int
  399. Width of one of 9 equal rectangles that will be used to compute
  400. a feature.
  401. height : int
  402. Height of one of 9 equal rectangles that will be used to compute
  403. a feature.
  404. lbp_code : int
  405. The descriptor of feature to visualize. If not provided, the
  406. descriptor with 0 value will be used.
  407. color_greater_block : tuple of 3 floats
  408. Floats specifying the color for the block that has greater
  409. intensity value. They should be in the range [0, 1].
  410. Corresponding values define (R, G, B) values. Default value
  411. is white (1, 1, 1).
  412. color_greater_block : tuple of 3 floats
  413. Floats specifying the color for the block that has greater intensity
  414. value. They should be in the range [0, 1]. Corresponding values define
  415. (R, G, B) values. Default value is cyan (0, 0.69, 0.96).
  416. alpha : float
  417. Value in the range [0, 1] that specifies opacity of visualization.
  418. 1 - fully transparent, 0 - opaque.
  419. Returns
  420. -------
  421. output : ndarray of float
  422. Image with MB-LBP visualization.
  423. References
  424. ----------
  425. .. [1] L. Zhang, R. Chu, S. Xiang, S. Liao, S.Z. Li. "Face Detection Based
  426. on Multi-Block LBP Representation", In Proceedings: Advances in
  427. Biometrics, International Conference, ICB 2007, Seoul, Korea.
  428. http://www.cbsr.ia.ac.cn/users/scliao/papers/Zhang-ICB07-MBLBP.pdf
  429. :DOI:`10.1007/978-3-540-74549-5_2`
  430. """
  431. # Default colors for regions.
  432. # White is for the blocks that are brighter.
  433. # Cyan is for the blocks that has less intensity.
  434. color_greater_block = np.asarray(color_greater_block, dtype=np.float64)
  435. color_less_block = np.asarray(color_less_block, dtype=np.float64)
  436. # Copy array to avoid the changes to the original one.
  437. output = np.copy(image)
  438. # As the visualization uses RGB color we need 3 bands.
  439. if len(image.shape) < 3:
  440. output = gray2rgb(image)
  441. # Colors are specified in floats.
  442. output = img_as_float(output)
  443. # Offsets of neighbor rectangles relative to central one.
  444. # It has order starting from top left and going clockwise.
  445. neighbor_rect_offsets = (
  446. (-1, -1),
  447. (-1, 0),
  448. (-1, 1),
  449. (0, 1),
  450. (1, 1),
  451. (1, 0),
  452. (1, -1),
  453. (0, -1),
  454. )
  455. # Pre-multiply the offsets with width and height.
  456. neighbor_rect_offsets = np.array(neighbor_rect_offsets)
  457. neighbor_rect_offsets[:, 0] *= height
  458. neighbor_rect_offsets[:, 1] *= width
  459. # Top-left coordinates of central rectangle.
  460. central_rect_r = r + height
  461. central_rect_c = c + width
  462. for element_num, offset in enumerate(neighbor_rect_offsets):
  463. offset_r, offset_c = offset
  464. curr_r = central_rect_r + offset_r
  465. curr_c = central_rect_c + offset_c
  466. has_greater_value = lbp_code & (1 << (7 - element_num))
  467. # Mix-in the visualization colors.
  468. if has_greater_value:
  469. new_value = (1 - alpha) * output[
  470. curr_r : curr_r + height, curr_c : curr_c + width
  471. ] + alpha * color_greater_block
  472. output[curr_r : curr_r + height, curr_c : curr_c + width] = new_value
  473. else:
  474. new_value = (1 - alpha) * output[
  475. curr_r : curr_r + height, curr_c : curr_c + width
  476. ] + alpha * color_less_block
  477. output[curr_r : curr_r + height, curr_c : curr_c + width] = new_value
  478. return output