diff --git a/doc/examples/plot_multiblock_local_binary_pattern.py b/doc/examples/plot_multiblock_local_binary_pattern.py new file mode 100644 index 00000000..f54d944e --- /dev/null +++ b/doc/examples/plot_multiblock_local_binary_pattern.py @@ -0,0 +1,73 @@ +""" +=========================================================== +Multi-Block Local Binary Pattern for texture classification +=========================================================== + +This example shows how to compute multi-block local binary pattern (MB-LBP) +features as well as how to visualize them. + +The features are calculated similarly to local binary patterns (LBPs), except +that summed blocks are used instead of individual pixel values. + +MB-LBP is an extension of LBP that can be computed on multiple scales in +constant time using the integral image. 9 equally-sized rectangles are used to +compute a feature. For each rectangle, the sum of the pixel intensities is +computed. Comparisons of these sums to that of the central rectangle determine +the feature, similarly to LBP (See `LBP `_). + +First, we generate an image to illustrate the functioning of MB-LBP: consider +a (9, 9) rectangle and divide it into (3, 3) block, upon which we then apply +MB-LBP. + +""" +from __future__ import print_function +from skimage.feature import multiblock_lbp +import numpy as np +from numpy.testing import assert_equal +from skimage.transform import integral_image + +# Create test matrix where first and fifth rectangles starting +# from top left clockwise have greater value than the central one. +test_img = np.zeros((9, 9), dtype='uint8') +test_img[3:6, 3:6] = 1 +test_img[:3, :3] = 50 +test_img[6:, 6:] = 50 + +# First and fifth bits should be filled. This correct value will +# be compared to the computed one. +correct_answer = 0b10001000 + +int_img = integral_image(test_img) + +lbp_code = multiblock_lbp(int_img, 0, 0, 3, 3) + +assert_equal(correct_answer, lbp_code) + +""" +Now let's apply the operator to a real image and see how the +visualization works. +""" +from skimage import data +from matplotlib import pyplot as plt +from skimage.feature import draw_multiblock_lbp + +test_img = data.coins() + +int_img = integral_image(test_img) + +lbp_code = multiblock_lbp(int_img, 0, 0, 90, 90) + +img = draw_multiblock_lbp(test_img, 0, 0, 90, 90, + lbp_code=lbp_code, alpha=0.5) + + +plt.imshow(img, interpolation='nearest') + +""" +.. image:: PLOT2RST.current_figure + +On the above plot we see the result of computing a MB-LBP and visualization of +the computed feature. The rectangles that have less intensities' sum than the +central rectangle are marked in cyan. The ones that have higher intensity +values are marked in white. The central rectangle is left untouched. +""" diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 567f561f..4d40e2c8 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1,7 +1,11 @@ from ._canny import canny from ._daisy import daisy from ._hog import hog -from .texture import greycomatrix, greycoprops, local_binary_pattern +from .texture import (greycomatrix, greycoprops, + local_binary_pattern, + multiblock_lbp, + draw_multiblock_lbp) + from .peak import peak_local_max from .corner import (corner_kitchen_rosenfeld, corner_harris, corner_shi_tomasi, corner_foerstner, corner_subpix, @@ -25,6 +29,8 @@ __all__ = ['canny', 'greycomatrix', 'greycoprops', 'local_binary_pattern', + 'multiblock_lbp', + 'draw_multiblock_lbp', 'peak_local_max', 'structure_tensor', 'structure_tensor_eigvals', diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx index d716e690..b90ab6dd 100644 --- a/skimage/feature/_texture.pyx +++ b/skimage/feature/_texture.pyx @@ -6,6 +6,7 @@ import numpy as np cimport numpy as cnp from libc.math cimport sin, cos, abs from .._shared.interpolation cimport bilinear_interpolation, round +from .._shared.transform cimport integrate cdef extern from "numpy/npy_math.h": @@ -264,3 +265,94 @@ def _local_binary_pattern(double[:, ::1] image, output[r, c] = lbp return np.asarray(output) + + +# Constant values that are used by `_multiblock_lbp` function. +# Values represent offsets of neighbour 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]) + Py_ssize_t[::1] mlbp_c_offsets = np.asarray([-1, 0, 1, 1, 1, 0, -1, -1]) + + +def _multiblock_lbp(float[:, ::1] int_image, + Py_ssize_t r, + Py_ssize_t c, + Py_ssize_t width, + Py_ssize_t height): + """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] Face Detection Based on Multi-Block LBP + Representation. Lun Zhang, Rufeng Chu, Shiming Xiang, Shengcai Liao, + Stan Z. Li + http://www.cbsr.ia.ac.cn/users/scliao/papers/Zhang-ICB07-MBLBP.pdf + """ + + 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 + + # Copy offset array to multiply it by width and height later. + Py_ssize_t[::1] r_offsets = mlbp_r_offsets.copy() + Py_ssize_t[::1] c_offsets = mlbp_c_offsets.copy() + + Py_ssize_t current_rect_r, current_rect_c + Py_ssize_t element_num, i + double current_rect_val + int has_greater_value + int lbp_code = 0 + + # Pre-multiply offsets with width and height. + for i in range(8): + r_offsets[i] = r_offsets[i]*height + c_offsets[i] = c_offsets[i]*width + + # Sum of intensity values of central rectangle. + cdef float 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 + r_offsets[element_num] + current_rect_c = central_rect_c + c_offsets[element_num] + + + 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 diff --git a/skimage/feature/tests/test_texture.py b/skimage/feature/tests/test_texture.py index 10444aa8..ae54a6bb 100644 --- a/skimage/feature/tests/test_texture.py +++ b/skimage/feature/tests/test_texture.py @@ -1,7 +1,11 @@ import numpy as np -from skimage.feature import greycomatrix, greycoprops, local_binary_pattern -from skimage._shared.testing import test_parallel +from skimage.feature import (greycomatrix, + greycoprops, + local_binary_pattern, + multiblock_lbp) +from skimage._shared.testing import test_parallel +from skimage.transform import integral_image class TestGLCM(): @@ -229,5 +233,28 @@ class TestLBP(): np.testing.assert_array_almost_equal(lbp, ref) +class TestMBLBP(): + + def test_single_mblbp(self): + + # Create dummy matrix where first and fifth rectangles have greater + # value than the central one. Therefore, the following bits + # should be 1. + test_img = np.zeros((9, 9), dtype='uint8') + test_img[3:6, 3:6] = 1 + test_img[:3, :3] = 255 + test_img[6:, 6:] = 255 + + # MB-LBP is filled in reverse order. So the first and fifth bits from + # the end should be filled. + correct_answer = 0b10001000 + + int_img = integral_image(test_img) + + lbp_code = multiblock_lbp(int_img, 0, 0, 3, 3) + + np.testing.assert_equal(lbp_code, correct_answer) + + if __name__ == '__main__': np.testing.run_module_suite() diff --git a/skimage/feature/texture.py b/skimage/feature/texture.py index beb66ff9..c26eb59d 100644 --- a/skimage/feature/texture.py +++ b/skimage/feature/texture.py @@ -4,7 +4,11 @@ Methods to characterize image textures. import numpy as np from .._shared.utils import assert_nD -from ._texture import _glcm_loop, _local_binary_pattern +from ..util import img_as_float +from ..color import gray2rgb +from ._texture import (_glcm_loop, + _local_binary_pattern, + _multiblock_lbp) def greycomatrix(image, distances, angles, levels=256, symmetric=False, @@ -291,3 +295,157 @@ def local_binary_pattern(image, P, R, method='default'): image = np.ascontiguousarray(image, dtype=np.double) output = _local_binary_pattern(image, P, R, methods[method.lower()]) return output + + +def multiblock_lbp(int_image, r, c, width, height): + """Multi-block local binary pattern (MB-LBP) [1]_. + + The features are calculated similarly to local binary patterns (LBPs), + (See :py:meth:`local_binary_pattern`) except that summed blocks are + used instead of individual pixel values. + + MB-LBP is an extension of LBP that can be computed on multiple scales + in constant time using the integral image. Nine equally-sized rectangles + are used to compute a feature. For each rectangle, the sum of the pixel + intensities is computed. Comparisons of these sums to that of the central + rectangle determine the feature, similarly to LBP. + + Parameters + ---------- + int_image : (N, M) 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 the 9 equal rectangles that will be used to compute + a feature. + height : int + Height of one of the 9 equal rectangles that will be used to compute + a feature. + + Returns + ------- + output : int + 8-bit MB-LBP feature descriptor. + + References + ---------- + .. [1] Face Detection Based on Multi-Block LBP + Representation. Lun Zhang, Rufeng Chu, Shiming Xiang, Shengcai Liao, + Stan Z. Li + http://www.cbsr.ia.ac.cn/users/scliao/papers/Zhang-ICB07-MBLBP.pdf + """ + + int_image = np.ascontiguousarray(int_image, dtype=np.float32) + lbp_code = _multiblock_lbp(int_image, r, c, width, height) + return lbp_code + + +def draw_multiblock_lbp(img, r, c, width, height, + lbp_code=0, + color_greater_block=[1, 1, 1], + color_less_block=[0, 0.69, 0.96], + alpha=0.5 + ): + """Multi-block local binary pattern visualization. + + Blocks with higher sums are colored with alpha-blended white rectangles, + whereas blocks with lower sums are colored alpha-blended cyan. Colors + and the `alpha` parameter can be changed. + + Parameters + ---------- + img : ndarray of float or uint + Image on which to visualize the pattern. + 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. + lbp_code : int + The descriptor of feature to visualize. If not provided, the + descriptor with 0 value will be used. + color_greater_block : list of 3 floats + Floats specifying the color for the block that has greater + intensity value. They should be in the range [0, 1]. + Corresponding values define (R, G, B) values. Default value + is white [1, 1, 1]. + color_greater_block : list of 3 floats + Floats specifying the color for the block that has greater intensity + value. They should be in the range [0, 1]. Corresponding values define + (R, G, B) values. Default value is cyan [0, 0.69, 0.96]. + alpha : float + Value in the range [0, 1] that specifies opacity of visualization. + 1 - fully transparent, 0 - opaque. + + Returns + ------- + output : ndarray of float + Image with MB-LBP visualization. + + References + ---------- + .. [1] Face Detection Based on Multi-Block LBP + Representation. Lun Zhang, Rufeng Chu, Shiming Xiang, Shengcai Liao, + Stan Z. Li + http://www.cbsr.ia.ac.cn/users/scliao/papers/Zhang-ICB07-MBLBP.pdf + """ + + # Default colors for regions. + # White is for the blocks that are brighter. + # Cyan is for the blocks that has less intensity. + color_greater_block = np.asarray(color_greater_block, dtype=np.float64) + color_less_block = np.asarray(color_less_block, dtype=np.float64) + + # Copy array to avoid the changes to the original one. + output = np.copy(img) + + # As the visualization uses RGB color we need 3 bands. + if len(img.shape) < 3: + output = gray2rgb(img) + + # Colors are specified in floats. + output = img_as_float(output) + + # Offsets of neighbour rectangles relative to central one. + # It has order starting from top left and going clockwise. + neighbour_rect_offsets = ((-1, -1), (-1, 0), (-1, 1), + (0, 1), (1, 1), (1, 0), + (1, -1), (0, -1)) + + # Pre-multiply the offsets with width and height. + neighbour_rect_offsets = np.array(neighbour_rect_offsets) + neighbour_rect_offsets[:, 0] *= height + neighbour_rect_offsets[:, 1] *= width + + # Top-left coordinates of central rectangle. + central_rect_r = r + height + central_rect_c = c + width + + for element_num, offset in enumerate(neighbour_rect_offsets): + + offset_r, offset_c = offset + + curr_r = central_rect_r + offset_r + curr_c = central_rect_c + offset_c + + has_greater_value = lbp_code & (1 << (7-element_num)) + + # Mix-in the visualization colors. + if has_greater_value: + new_value = ((1-alpha) * output[curr_r:curr_r+height, curr_c:curr_c+width] + + alpha * color_greater_block) + output[curr_r:curr_r+height, curr_c:curr_c+width] = new_value + else: + new_value = ((1-alpha) * output[curr_r:curr_r+height, curr_c:curr_c+width] + + alpha * color_less_block) + output[curr_r:curr_r+height, curr_c:curr_c+width] = new_value + + return output