diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 735bffa9..e909a0a0 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -108,6 +108,7 @@ Adaptive thresholding Implementation of Matlab's `regionprops` Estimation of geometric transformation parameters - + Local binary pattern texture classification + - Pavel Campr Fixes and tests for Histograms of Oriented Gradients. diff --git a/doc/examples/plot_local_binary_pattern.py b/doc/examples/plot_local_binary_pattern.py new file mode 100644 index 00000000..85fd5b95 --- /dev/null +++ b/doc/examples/plot_local_binary_pattern.py @@ -0,0 +1,87 @@ +""" +=============================================== +Local Binary Pattern for texture classification +=============================================== + +In this example, we will see how to classify textures based on LBP (Local +Binary Pattern). The histogram of the LBP result is a good measure to classify +textures. For simplicity the histogram distributions are then tested against +each other using the Kullback-Leibler-Divergence. +""" + +import numpy as np +import matplotlib +import matplotlib.pyplot as plt +import scipy.ndimage as nd +import skimage.feature as ft +from skimage import data + + +# settings for LBP +METHOD = 'uniform' +P = 16 +R = 2 +matplotlib.rcParams['font.size'] = 9 + + +def kullback_leibler_divergence(p, q): + p = np.asarray(p) + q = np.asarray(q) + filt = np.logical_and(p != 0, q != 0) + return np.sum(p[filt] * np.log2(p[filt] / q[filt])) + + +def match(refs, img): + best_score = 10 + best_name = None + lbp = ft.local_binary_pattern(img, P, R, METHOD) + hist, _ = np.histogram(lbp, normed=True, bins=P + 2, range=(0, P + 2)) + for name, ref in refs.items(): + ref_hist, _ = np.histogram(ref, normed=True, bins=P + 2, + range=(0, P + 2)) + score = kullback_leibler_divergence(hist, ref_hist) + if score < best_score: + best_score = score + best_name = name + return best_name + + +brick = data.load('brick.png') +grass = data.load('grass.png') +wall = data.load('rough-wall.png') + +refs = { + 'brick': ft.local_binary_pattern(brick, P, R, METHOD), + 'grass': ft.local_binary_pattern(grass, P, R, METHOD), + 'wall': ft.local_binary_pattern(wall, P, R, METHOD) +} + +# classify rotated textures +print 'Rotated images matched against references using LBP:' +print 'original: brick, rotated: 30deg, match result:', +print match(refs, nd.rotate(brick, angle=30, reshape=False)) +print 'original: brick, rotated: 70deg, match result:', +print match(refs, nd.rotate(brick, angle=70, reshape=False)) +print 'original: grass, rotated: 145deg, match result:', +print match(refs, nd.rotate(grass, angle=145, reshape=False)) + +# plot histograms of LBP of textures +fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows=2, ncols=3, + figsize=(9, 6)) +plt.gray() + +ax1.imshow(brick) +ax1.axis('off') +ax4.hist(refs['brick'].ravel(), normed=True, bins=P + 2, range=(0, P + 2)) +ax4.set_ylabel('Percentage') + +ax2.imshow(grass) +ax2.axis('off') +ax5.hist(refs['grass'].ravel(), normed=True, bins=P + 2, range=(0, P + 2)) +ax5.set_xlabel('Uniform LBP values') + +ax3.imshow(wall) +ax3.axis('off') +ax6.hist(refs['wall'].ravel(), normed=True, bins=P + 2, range=(0, P + 2)) + +plt.show() diff --git a/skimage/data/brick.png b/skimage/data/brick.png new file mode 100644 index 00000000..c69e71b2 Binary files /dev/null and b/skimage/data/brick.png differ diff --git a/skimage/data/grass.png b/skimage/data/grass.png new file mode 100644 index 00000000..16062a35 Binary files /dev/null and b/skimage/data/grass.png differ diff --git a/skimage/data/rough-wall.png b/skimage/data/rough-wall.png new file mode 100644 index 00000000..8e2ced60 Binary files /dev/null and b/skimage/data/rough-wall.png differ diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index eb98945c..1597e9a8 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1,5 +1,5 @@ from ._hog import hog -from ._greycomatrix import greycomatrix, greycoprops +from .texture import greycomatrix, greycoprops, local_binary_pattern from .peak import peak_local_max from ._harris import harris from .template import match_template diff --git a/skimage/feature/_greycomatrix_cy.pyx b/skimage/feature/_greycomatrix_cy.pyx deleted file mode 100644 index b9ff7f23..00000000 --- a/skimage/feature/_greycomatrix_cy.pyx +++ /dev/null @@ -1,67 +0,0 @@ -"""Cython implementation for computing a grey level co-occurance matrix -""" - -import numpy as np -cimport numpy as np -cimport cython - -cdef extern from "math.h": - double sin(double) - double cos(double) - -@cython.boundscheck(False) -def _glcm_loop(np.ndarray[dtype=np.uint8_t, ndim=2, - negative_indices=False, mode='c'] image, - np.ndarray[dtype=np.float64_t, ndim=1, - negative_indices=False, mode='c'] distances, - np.ndarray[dtype=np.float64_t, ndim=1, - negative_indices=False, mode='c'] angles, - int levels, - np.ndarray[dtype=np.uint32_t, ndim=4, - negative_indices=False, mode='c'] out - ): - """Perform co-occurnace matrix accumulation - - Parameters - ---------- - image : ndarray - Input image, which is converted to the uint8 data type. - 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 grey-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: - np.int32_t a_inx, d_idx - np.int32_t r, c, rows, cols, row, col - np.int32_t i, j - - rows = image.shape[0] - cols = image.shape[1] - - for a_idx, angle in enumerate(angles): - for d_idx, distance in enumerate(distances): - for r in range(rows): - for c in range(cols): - i = image[r, c] - - # compute the location of the offset pixel - row = r + (sin(angle) * distance + 0.5) - col = c + (cos(angle) * distance + 0.5); - - # make sure the offset is within bounds - if row >= 0 and row < rows and \ - col >= 0 and col < cols: - j = image[row, col] - - if i >= 0 and i < levels and \ - j >= 0 and j < levels: - out[i, j, d_idx, a_idx] += 1 diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx new file mode 100644 index 00000000..ca51e177 --- /dev/null +++ b/skimage/feature/_texture.pyx @@ -0,0 +1,181 @@ +#cython: cdivison=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False +import numpy as np +cimport numpy as np +from libc.math cimport sin, cos, abs +from skimage.transform._project cimport bilinear_interpolation + + +def _glcm_loop(np.ndarray[dtype=np.uint8_t, ndim=2, + negative_indices=False, mode='c'] image, + np.ndarray[dtype=np.float64_t, ndim=1, + negative_indices=False, mode='c'] distances, + np.ndarray[dtype=np.float64_t, ndim=1, + negative_indices=False, mode='c'] angles, + int levels, + np.ndarray[dtype=np.uint32_t, ndim=4, + negative_indices=False, mode='c'] out + ): + """Perform co-occurrence matrix accumulation. + + Parameters + ---------- + image : ndarray + Input image, which is converted to the uint8 data type. + 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 grey-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: + np.int32_t a_inx, d_idx + np.int32_t r, c, rows, cols, row, col + np.int32_t i, j + + rows = image.shape[0] + cols = image.shape[1] + + for a_idx, angle in enumerate(angles): + for d_idx, distance in enumerate(distances): + for r in range(rows): + for c in range(cols): + i = image[r, c] + + # compute the location of the offset pixel + row = r + (sin(angle) * distance + 0.5) + col = c + (cos(angle) * distance + 0.5); + + # make sure the offset is within bounds + if row >= 0 and row < rows and \ + col >= 0 and col < cols: + j = image[row, col] + + if i >= 0 and i < levels and \ + j >= 0 and j < levels: + out[i, j, d_idx, a_idx] += 1 + + +cdef inline int _bit_rotate_right(int value, int length): + """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(np.ndarray[double, ndim=2] image, + int P, float R, char method='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) double array + Graylevel image. + P : int + Number of circularly symmetric neighbour set points (quantization of the + angular space). + R : float + Radius of circle (spatial resolution of the operator). + method : {'D', 'R', 'U', 'V'} + Method to determine the pattern:: + * 'D': 'default' + * 'R': 'ror' + * 'U': 'uniform' + * 'V': 'var' + + Returns + ------- + output : (N, M) array + LBP image. + """ + + # texture weights + cdef np.ndarray[int, ndim=1] weights = 2 ** np.arange(P, dtype=np.int32) + # local position of texture elements + rp = - R * np.sin(2 * np.pi * np.arange(P, dtype=np.double) / P) + cp = R * np.cos(2 * np.pi * np.arange(P, dtype=np.double) / P) + cdef np.ndarray[double, ndim=2] coords = np.round(np.vstack([rp, cp]).T, 5) + + # pre allocate arrays for computation + cdef np.ndarray[double, ndim=1] texture = np.zeros(P, np.double) + cdef np.ndarray[char, ndim=1] signed_texture = np.zeros(P, np.int8) + cdef np.ndarray[int, ndim=1] rotation_chain = np.zeros(P, np.int32) + + output_shape = (image.shape[0], image.shape[1]) + cdef np.ndarray[double, ndim=2] output = np.zeros(output_shape, np.double) + + cdef int rows = image.shape[0] + cdef int cols = image.shape[1] + + cdef double lbp + cdef int r, c, changes, i + for r in range(image.shape[0]): + for c in range(image.shape[1]): + for i in range(P): + texture[i] = bilinear_interpolation(image.data, + rows, cols, r + coords[i, 0], c + coords[i, 1], 'C') + # 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 == 'uniform' or method == 'var': + if method == 'U' or method == 'V': + # determine number of 0 - 1 changes + changes = 0 + for i in range(P - 1): + changes += abs(signed_texture[i] - signed_texture[i + 1]) + + if changes <= 2: + for i in range(P): + lbp += signed_texture[i] + else: + lbp = P + 1 + + if method == 'V': + var = np.var(texture) + if var != 0: + lbp /= var + else: + lbp = np.nan + else: + # method == 'default' + for i in range(P): + lbp += signed_texture[i] * weights[i] + + # method == 'ror' + if method == 'R': + # shift LBP P times to the right and get minimum value + rotation_chain[0] = 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 output diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py index 2c50a592..0b0b80bd 100644 --- a/skimage/feature/setup.py +++ b/skimage/feature/setup.py @@ -12,11 +12,12 @@ def configuration(parent_package='', top_path=None): config = Configuration('feature', parent_package, top_path) config.add_data_dir('tests') - cython(['_greycomatrix_cy.pyx'], working_path=base_path) + cython(['_texture.pyx'], working_path=base_path) cython(['_template.pyx'], working_path=base_path) - config.add_extension('_greycomatrix_cy', sources=['_greycomatrix_cy.c'], - include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_texture', sources=['_texture.c'], + include_dirs=[get_numpy_include_dirs(), + '../transform']) config.add_extension('_template', sources=['_template.c'], include_dirs=[get_numpy_include_dirs()]) diff --git a/skimage/feature/tests/test_glcm.py b/skimage/feature/tests/test_texture.py similarity index 69% rename from skimage/feature/tests/test_glcm.py rename to skimage/feature/tests/test_texture.py index da90ec56..d48a14f7 100644 --- a/skimage/feature/tests/test_glcm.py +++ b/skimage/feature/tests/test_texture.py @@ -1,8 +1,9 @@ import numpy as np -from skimage.feature import greycomatrix, greycoprops +from skimage.feature import greycomatrix, greycoprops, local_binary_pattern class TestGLCM(): + def setup(self): self.image = np.array([[0, 0, 1, 1], [0, 0, 1, 1], @@ -140,5 +141,63 @@ class TestGLCM(): 'energy', 'correlation', 'ASM']: greycoprops(result, prop) + +class TestLBP(): + + def setup(self): + self.image = np.array([[255, 6, 255, 0, 141, 0], + [ 48, 250, 204, 166, 223, 63], + [ 8, 0, 159, 50, 255, 30], + [167, 255, 63, 40, 128, 255], + [ 0, 255, 30, 34, 255, 24], + [146, 241, 255, 0, 189, 126]], dtype='double') + + def test_default(self): + lbp = local_binary_pattern(self.image, 8, 1, 'default') + ref = np.array([[ 0, 251, 0, 255, 96, 255], + [143, 0, 20, 153, 64, 56], + [238, 255, 12, 191, 0, 252], + [129, 64., 62, 159, 199, 0], + [255, 4, 255, 175, 0, 254], + [ 3, 5, 0, 255, 4, 24]]) + np.testing.assert_array_equal(lbp, ref) + + def test_ror(self): + lbp = local_binary_pattern(self.image, 8, 1, 'ror') + ref = np.array([[ 0, 127, 0, 255, 3, 255], + [ 31, 0, 5, 51, 1, 7], + [119, 255, 3, 127, 0, 63], + [ 3, 1, 31, 63, 31, 0], + [255, 1, 255, 95, 0, 127], + [ 3, 5, 0, 255, 1, 3]]) + np.testing.assert_array_equal(lbp, ref) + + def test_uniform(self): + lbp = local_binary_pattern(self.image, 8, 1, 'uniform') + ref = np.array([[0, 7, 0, 8, 2, 8], + [5, 0, 9, 9, 1, 3], + [9, 8, 2, 7, 0, 6], + [2, 1, 5, 6, 5, 0], + [8, 1, 8, 9, 0, 7], + [2, 9, 0, 8, 1, 2]]) + np.testing.assert_array_equal(lbp, ref) + + def test_var(self): + lbp = local_binary_pattern(self.image, 8, 1, 'var') + ref = np.array([[0. , 0.00072786, 0. , 0.00115377, + 0.00032355, 0.00224467], + [0.00051758, 0. , 0.0026383 , 0.00163246, + 0.00027414, 0.00041124], + [0.00192834, 0.00130368, 0.00042095, 0.00171894, + 0. , 0.00063726], + [0.00023048, 0.00019464 , 0.00082291, 0.00225386, + 0.00076696, 0. ], + [0.00097253, 0.00013236, 0.0009134 , 0.0014467 , + 0. , 0.00082472], + [0.00024701, 0.0012277 , 0. , 0.00109869, + 0.00015445, 0.00035881]]) + np.testing.assert_array_almost_equal(lbp, ref) + + if __name__ == '__main__': np.testing.run_module_suite() diff --git a/skimage/feature/texture.py b/skimage/feature/texture.py new file mode 100644 index 00000000..aa970d18 --- /dev/null +++ b/skimage/feature/texture.py @@ -0,0 +1,278 @@ +""" +Methods to characterize image textures. +""" + +import math +import numpy as np +from scipy import ndimage + +from ._texture import _glcm_loop, _local_binary_pattern + + +def greycomatrix(image, distances, angles, levels=256, symmetric=False, + normed=False): + """Calculate the grey-level co-occurrence matrix. + + A grey level co-occurence matrix is a histogram of co-occuring + greyscale values at a given offset over an image. + + Parameters + ---------- + image : array_like of uint8 + Integer typed input image. The image will be cast to uint8, so + the maximum value must be less than 256. + distances : array_like + List of pixel pair distance offsets. + angles : array_like + List of pixel pair angles in radians. + levels : int, optional + The input image should contain integers in [0, levels-1], + where levels indicate the number of grey-levels counted + (typically 256 for an 8-bit image). The maximum value is + 256. + symmetric : bool, optional + If True, the output matrix `P[:, :, d, theta]` is symmetric. This + is accomplished by ignoring the order of value pairs, so both + (i, j) and (j, i) are accumulated when (i, j) is encountered + for a given offset. The default is False. + normed : bool, optional + If True, normalize each matrix `P[:, :, d, theta]` by dividing + by the total number of accumulated co-occurrences for the given + offset. The elements of the resulting matrix sum to 1. The + default is False. + + Returns + ------- + P : 4-D ndarray + The grey-level co-occurrence histogram. The value + `P[i,j,d,theta]` is the number of times that grey-level `j` + occurs at a distance `d` and at an angle `theta` from + grey-level `i`. If `normed` is `False`, the output is of + type uint32, otherwise it is float64. + + References + ---------- + .. [1] The GLCM Tutorial Home Page, + http://www.fp.ucalgary.ca/mhallbey/tutorial.htm + .. [2] Pattern Recognition Engineering, Morton Nadler & Eric P. + Smith + .. [3] Wikipedia, http://en.wikipedia.org/wiki/Co-occurrence_matrix + + + Examples + -------- + Compute 2 GLCMs: One for a 1-pixel offset to the right, and one + for a 1-pixel offset upwards. + + >>> image = np.array([[0, 0, 1, 1], + ... [0, 0, 1, 1], + ... [0, 2, 2, 2], + ... [2, 2, 3, 3]], dtype=np.uint8) + >>> result = greycomatrix(image, [1], [0, np.pi/2], levels=4) + >>> result[:, :, 0, 0] + array([[2, 2, 1, 0], + [0, 2, 0, 0], + [0, 0, 3, 1], + [0, 0, 0, 1]], dtype=uint32) + >>> result[:, :, 0, 1] + array([[3, 0, 2, 0], + [0, 2, 2, 0], + [0, 0, 1, 2], + [0, 0, 0, 0]], dtype=uint32) + + """ + + assert levels <= 256 + image = np.ascontiguousarray(image) + assert image.ndim == 2 + assert image.min() >= 0 + assert image.max() < levels + image = image.astype(np.uint8) + distances = np.ascontiguousarray(distances, dtype=np.float64) + angles = np.ascontiguousarray(angles, dtype=np.float64) + assert distances.ndim == 1 + assert angles.ndim == 1 + + P = np.zeros((levels, levels, len(distances), len(angles)), + dtype=np.uint32, order='C') + + # count co-occurences + _glcm_loop(image, distances, angles, levels, P) + + # make each GLMC symmetric + if symmetric: + Pt = np.transpose(P, (1, 0, 2, 3)) + P = P + Pt + + # normalize each GLMC + if normed: + P = P.astype(np.float64) + glcm_sums = np.apply_over_axes(np.sum, P, axes=(0, 1)) + glcm_sums[glcm_sums == 0] = 1 + P /= glcm_sums + + return P + + +def greycoprops(P, prop='contrast'): + """Calculate texture properties of a GLCM. + + Compute a feature of a grey level co-occurrence matrix to serve as + a compact summary of the matrix. The properties are computed as + follows: + + - 'contrast': :math:`\\sum_{i,j=0}^{levels-1} P_{i,j}(i-j)^2` + - 'dissimilarity': :math:`\\sum_{i,j=0}^{levels-1}P_{i,j}|i-j|` + - 'homogeneity': :math:`\\sum_{i,j=0}^{levels-1}\\frac{P_{i,j}}{1+(i-j)^2}` + - 'ASM': :math:`\\sum_{i,j=0}^{levels-1} P_{i,j}^2` + - 'energy': :math:`\\sqrt{ASM}` + - 'correlation': + .. math:: \\sum_{i,j=0}^{levels-1} P_{i,j}\\left[\\frac{(i-\\mu_i) \\ + (j-\\mu_j)}{\\sqrt{(\\sigma_i^2)(\\sigma_j^2)}}\\right] + + + Parameters + ---------- + P : ndarray + Input array. `P` is the grey-level co-occurrence histogram + for which to compute the specified property. The value + `P[i,j,d,theta]` is the number of times that grey-level j + occurs at a distance d and at an angle theta from + grey-level i. + + prop : {'contrast', 'dissimilarity', 'homogeneity', 'energy', \ + 'correlation', 'ASM'}, optional + The property of the GLCM to compute. The default is 'contrast'. + + Returns + ------- + results : 2-D ndarray + 2-dimensional array. `results[d, a]` is the property 'prop' for + the d'th distance and the a'th angle. + + References + ---------- + .. [1] The GLCM Tutorial Home Page, + http://www.fp.ucalgary.ca/mhallbey/tutorial.htm + + Examples + -------- + Compute the contrast for GLCMs with distances [1, 2] and angles + [0 degrees, 90 degrees] + + >>> image = np.array([[0, 0, 1, 1], + ... [0, 0, 1, 1], + ... [0, 2, 2, 2], + ... [2, 2, 3, 3]], dtype=np.uint8) + >>> g = greycomatrix(image, [1, 2], [0, np.pi/2], levels=4, + ... normed=True, symmetric=True) + >>> contrast = greycoprops(g, 'contrast') + >>> contrast + array([[ 0.58333333, 1. ], + [ 1.25 , 2.75 ]]) + + """ + + assert P.ndim == 4 + (num_level, num_level2, num_dist, num_angle) = P.shape + assert num_level == num_level2 + assert num_dist > 0 + assert num_angle > 0 + + # create weights for specified property + I, J = np.ogrid[0:num_level, 0:num_level] + if prop == 'contrast': + weights = (I - J) ** 2 + elif prop == 'dissimilarity': + weights = np.abs(I - J) + elif prop == 'homogeneity': + weights = 1. / (1. + (I - J) ** 2) + elif prop in ['ASM', 'energy', 'correlation']: + pass + else: + raise ValueError('%s is an invalid property' % (prop)) + + # compute property for each GLCM + if prop == 'energy': + asm = np.apply_over_axes(np.sum, (P ** 2), axes=(0, 1))[0, 0] + results = np.sqrt(asm) + elif prop == 'ASM': + results = np.apply_over_axes(np.sum, (P ** 2), axes=(0, 1))[0, 0] + elif prop == 'correlation': + results = np.zeros((num_dist, num_angle), dtype=np.float64) + I = np.array(range(num_level)).reshape((num_level, 1, 1, 1)) + J = np.array(range(num_level)).reshape((1, num_level, 1, 1)) + diff_i = I - np.apply_over_axes(np.sum, (I * P), axes=(0, 1))[0, 0] + diff_j = J - np.apply_over_axes(np.sum, (J * P), axes=(0, 1))[0, 0] + + std_i = np.sqrt(np.apply_over_axes(np.sum, (P * (diff_i) ** 2), + axes=(0, 1))[0, 0]) + std_j = np.sqrt(np.apply_over_axes(np.sum, (P * (diff_j) ** 2), + axes=(0, 1))[0, 0]) + cov = np.apply_over_axes(np.sum, (P * (diff_i * diff_j)), + axes=(0, 1))[0, 0] + + # handle the special case of standard deviations near zero + mask_0 = std_i < 1e-15 + mask_0[std_j < 1e-15] = True + results[mask_0] = 1 + + # handle the standard case + mask_1 = mask_0 == False + results[mask_1] = cov[mask_1] / (std_i[mask_1] * std_j[mask_1]) + elif prop in ['contrast', 'dissimilarity', 'homogeneity']: + weights = weights.reshape((num_level, num_level, 1, 1)) + results = np.apply_over_axes(np.sum, (P * weights), axes=(0, 1))[0, 0] + + return results + + +def local_binary_pattern(image, P, R, method='default'): + """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) array + Graylevel image. + P : int + Number of circularly symmetric neighbour set points (quantization of the + angular space). + R : float + Radius of circle (spatial resolution of the operator). + method : {'D', 'R', 'U', 'V'} + Method to determine the pattern:: + * 'default': original local binary pattern which is gray scale but not + rotation invariant. + * 'ror': extension of default implementation which is gray scale and + rotation invariant. + * 'uniform': improved rotation invariance with uniform patterns and + finer quantization of the angular space which is gray scale and + rotation invariant. + * 'var': rotation invariant variance measures of the contrast of local + image texture which is rotation but not gray scale invariant. + + Returns + ------- + output : (N, M) array + LBP image. + + References + ---------- + .. [1] Multiresolution Gray-Scale and Rotation Invariant Texture + Classification with Local Binary Patterns. + Timo Ojala, Matti Pietikainen, Topi Maenpaa. + http://www.rafbis.it/biplab15/images/stories/docenti/Danielriccio/\ + Articoliriferimento/LBP.pdf, 2002. + """ + + methods = { + 'default': ord('D'), + 'ror': ord('R'), + 'uniform': ord('U'), + 'var': ord('V') + } + image = np.array(image, dtype='double', copy=True) + output = _local_binary_pattern(image, P, R, methods[method.lower()]) + return output diff --git a/skimage/transform/_project.pxd b/skimage/transform/_project.pxd new file mode 100644 index 00000000..dee51026 --- /dev/null +++ b/skimage/transform/_project.pxd @@ -0,0 +1,7 @@ +cimport numpy as np +import numpy as np + + +cdef inline double bilinear_interpolation(double* image, int rows, int cols, + double r, double c, char mode, + double cval=*) diff --git a/skimage/transform/_project.pyx b/skimage/transform/_project.pyx index 734d6bff..3f31f229 100644 --- a/skimage/transform/_project.pyx +++ b/skimage/transform/_project.pyx @@ -1,38 +1,66 @@ -#cython: cdivison=True boundscheck=False +#cython: cdivison=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False -__all__ = ['homography'] - -cimport cython cimport numpy as np - import numpy as np -import cython - from cython.operator import dereference +from libc.math cimport ceil, floor -np.import_array() -cdef extern from "math.h": - double floor(double) - double fmod(double, double) +cdef inline double bilinear_interpolation(double* image, int rows, int cols, + double r, double c, char mode, + double cval=0): + """Bilinear interpolation at a given position in the image. -cdef double get_pixel(double *image, int rows, int cols, - int r, int c, char mode, double cval=0): + Parameters + ---------- + image : double array + Input image. + rows, cols: int + Shape of image. + r, c : int + Position at which to interpolate. + mode : {'C', 'W', 'M'} + Wrapping mode. Constant, Wrap or Mirror. + cval : double + Constant value to use for constant mode. + + """ + cdef double dr, dc + cdef int minr, minc, maxr, maxc + + minr = floor(r) + minc = floor(c) + maxr = ceil(r) + maxc = ceil(c) + dr = r - minr + dc = c - minc + top = (1 - dc) * get_pixel(image, rows, cols, minr, minc, mode, cval) \ + + dc * get_pixel(image, rows, cols, minr, maxc, mode, cval) + bottom = (1 - dc) * get_pixel(image, rows, cols, maxr, minc, mode, cval) \ + + dc * get_pixel(image, rows, cols, maxr, maxc, mode, cval) + return (1 - dr) * top + dr * bottom + + +cdef inline double get_pixel(double* image, int rows, int cols, int r, int c, + char mode, double cval=0): """Get a pixel from the image, taking wrapping mode into consideration. Parameters ---------- - image : *double + image : double array Input image. - rows, cols : int - Dimensions of image. + rows, cols: int + Shape of image. r, c : int Position at which to get the pixel. mode : {'C', 'W', 'M'} - Wrapping mode. Constant, Wrap or Mirror. + Wrapping mode. Constant, Wrap or Mirror. cval : double - Constant value to use for mode constant. - + Constant value to use for constant mode. + """ if mode == 'C': if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1): @@ -40,13 +68,13 @@ cdef double get_pixel(double *image, int rows, int cols, else: return image[r * cols + c] else: - return image[coord_map(rows, r, mode) * cols + - coord_map(cols, c, mode)] + return image[coord_map(rows, r, mode) * cols + coord_map(cols, c, mode)] -cdef int coord_map(int dim, int coord, char mode): + +cdef inline int coord_map(int dim, int coord, char mode): """ - Wrap a coordinate, according to a given dimension and mode. - + Wrap a coordinate, according to a given mode. + Parameters ---------- dim : int @@ -56,7 +84,7 @@ cdef int coord_map(int dim, int coord, char mode): mode : {'W', 'M'} Whether to wrap or mirror the coordinate if it falls outside [0, dim). - + """ dim = dim - 1 if mode == 'M': # mirror @@ -79,7 +107,9 @@ cdef int coord_map(int dim, int coord, char mode): return coord -cdef tf(double x, double y, double* H, double *x_, double *y_): + +cdef inline _matrix_transform(double x, double y, double* H, double *x_, + double *y_): """Apply a homography to a coordinate. Parameters @@ -98,18 +128,15 @@ cdef tf(double x, double y, double* H, double *x_, double *y_): yy = H[3] * x + H[4] * y + H[5] zz = H[6] * x + H[7] * y + H[8] - xx = xx / zz - yy = yy / zz + x_[0] = xx / zz + y_[0] = yy / zz - x_[0] = xx - y_[0] = yy -@cython.boundscheck(False) def homography(np.ndarray image, np.ndarray H, output_shape=None, mode='constant', double cval=0): """ Projective transformation (homography). - + Perform a projective transformation (homography) of a floating point image, using bi-linear interpolation. @@ -140,8 +167,6 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, Transformation matrix H that defines the homography. output_shape : tuple (rows, cols) Shape of the output image generated. - order : int - Order of splines used in interpolation. mode : {'constant', 'mirror', 'wrap'} How to handle values outside the image borders. cval : string @@ -150,8 +175,7 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, """ - cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] img = \ - np.ascontiguousarray(image, dtype=np.double) + cdef np.ndarray[dtype=np.double_t, ndim=2] img = image.astype(np.double) cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] M = \ np.ascontiguousarray(np.linalg.inv(H)) @@ -165,7 +189,6 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, elif mode == 'mirror': mode_c = ord('M') - cdef int out_r, out_c, columns, rows if output_shape is None: out_r = img.shape[0] out_c = img.shape[1] @@ -173,37 +196,18 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, out_r = output_shape[0] out_c = output_shape[1] - rows = img.shape[0] - columns = img.shape[1] - cdef np.ndarray[dtype=np.double_t, ndim=2] out = \ np.zeros((out_r, out_c), dtype=np.double) - - cdef int tfr, tfc, r_int, c_int - cdef double y0, y1, y2, y3 - cdef double r, c, z, t, u + + cdef int tfr, tfc + cdef double r, c + cdef int rows = img.shape[0] + cdef int cols = img.shape[1] for tfr in range(out_r): for tfc in range(out_c): - tf(tfc, tfr, M.data, &c, &r) - r_int = floor(r) - c_int = floor(c) - - t = r - r_int - u = c - c_int - - y0 = get_pixel(img.data, rows, columns, - r_int, c_int, mode_c) - y1 = get_pixel(img.data, rows, columns, - r_int + 1, c_int, mode_c) - y2 = get_pixel(img.data, rows, columns, - r_int + 1, c_int + 1, mode_c) - y3 = get_pixel(img.data, rows, columns, - r_int, c_int + 1, mode_c) - - out[tfr, tfc] = \ - (1 - t) * (1 - u) * y0 + \ - t * (1 - u) * y1 + \ - t * u * y2 + (1 - t) * u * y3; + _matrix_transform(tfc, tfr, M.data, &c, &r) + out[tfr, tfc] = bilinear_interpolation(img.data, rows, + cols, r, c, mode_c) return out diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py index 7c2e52f2..f35b7f57 100644 --- a/skimage/transform/tests/test_warps.py +++ b/skimage/transform/tests/test_warps.py @@ -2,7 +2,7 @@ from numpy.testing import assert_array_almost_equal, run_module_suite import numpy as np from skimage.transform import (warp, homography, fast_homography, - SimilarityTransform) + SimilarityTransform, ProjectiveTransform) from skimage import transform as tf, data, img_as_float from skimage.color import rgb2gray @@ -34,7 +34,7 @@ def test_homography(): def test_fast_homography(): - img = rgb2gray(data.lena()).astype(np.uint8) + img = rgb2gray(data.lena()) img = img[:, :100] theta = np.deg2rad(30) @@ -49,20 +49,19 @@ def test_fast_homography(): H[:2, 2] = [tx, ty] for mode in ('constant', 'mirror', 'wrap'): - p0 = homography(img, H, mode=mode, order=1) + p0 = warp(img, ProjectiveTransform(H).inverse, mode=mode, order=1) p1 = fast_homography(img, H, mode=mode) - p1 = np.round(p1) - ## import matplotlib.pyplot as plt - ## f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4) - ## ax0.imshow(img) - ## ax1.imshow(p0, cmap=plt.cm.gray) - ## ax2.imshow(p1, cmap=plt.cm.gray) - ## ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray) - ## plt.show() + # import matplotlib.pyplot as plt + # f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4) + # ax0.imshow(img) + # ax1.imshow(p0, cmap=plt.cm.gray) + # ax2.imshow(p1, cmap=plt.cm.gray) + # ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray) + # plt.show() d = np.mean(np.abs(p0 - p1)) - assert d < 0.2 + assert d < 0.001 def test_swirl():