wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-16 11:21:25 +08:00
Code Issues Packages Projects Releases Wiki Activity

add local binary pattern texture analysis

Browse Source
This commit is contained in:
Johannes Schönberger
2012-08-20 22:46:57 +02:00
parent 6410905b3b
commit d9444e7612
3 changed files with 406 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/feature/__init__.py
+2
View File
@@ -1,5 +1,7 @@
from ._hog import hog
from ._greycomatrix import greycomatrix, greycoprops
from .hog import hog
from ._texture import greycomatrix, greycoprops, local_binary_pattern
from .peak import peak_local_max
from ._harris import harris
from .template import match_template
skimage/feature/_texture.py
+340
View File
@@ -0,0 +1,340 @@
"""
Methods to characterize image textures.
"""
import math
import numpy as np
from scipy import ndimage
from ._greycomatrix import _glcm_loop
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 bit_rotate_right(value, 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(image, P, R, method='default'):
"""Texture classification using gray scale and rotation invariant LBP (Local
Binary Patterns).
Parameters
----------
image : NxM 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 : {'default', 'ror', 'uniform', 'var'}
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 : NxM array
LBP image
References
----------
Timo Ojala, Matti Pietikainen, Topi Maenpaa. Multiresolution Gray-Scale and
Rotation Invariant Texture Classification with Local Binary Patterns.
http://www.rafbis.it/biplab15/images/stories/docenti/Danielriccio/\
Articoliriferimento/LBP.pdf, 2002.
"""
method = method.lower()
# texture weights
weights = 2 ** np.arange(P)
# local position of texture elements
rp = - R * np.sin(2 * math.pi * np.arange(P) / P)
cp = R * np.cos(2 * math.pi * np.arange(P) / P)
coords = np.vstack([rp, cp]) + math.ceil(R)
# maximum size of neighbourhood for filtering
max_size = 2 * math.ceil(R) + 1
# center index of flattened neightbourhood
center_index = (max_size ** 2 - 1) / 2
if method == 'ror':
# allocate array for rotation invariance
rotation_chain = np.zeros(P, dtype='int')
def compute_lbp(texture):
# subtract value of center pixel
texture -= texture[center_index]
#: get texture elements using bilinear interpolation
texture = texture.reshape(max_size, max_size)
texture = ndimage.map_coordinates(texture, coords, order=1)
#: signed / thresholded texture
signed = texture.copy()
signed[signed>=0] = 1
signed[signed<0] = 0
if method in ('uniform', 'var'):
#: determine number of 0 - 1 changes
changes = np.sum(np.abs(np.diff(signed)))
if changes <= 2:
lbp = np.sum(signed)
else:
lbp = P + 1
if method == 'var':
lbp /= np.var(texture)
else:
# method == 'default'
lbp = np.sum(signed * weights)
if method == 'ror':
#: shift LBP P times to the right and get minimum value
rotation_chain[0] = lbp
for i in xrange(1, P):
rotation_chain[i] = bit_rotate_right(rotation_chain[i-1], P)
lbp = np.min(rotation_chain)
return lbp
dtype = 'int'
if method == 'var':
dtype = 'float'
output = np.zeros(image.shape, dtype)
ndimage.generic_filter(image, compute_lbp, size=(max_size, max_size),
mode='constant', cval=0, output=output)
return output
skimage/feature/tests/test_glcm.py → skimage/feature/tests/test_texture.py
+64 -1
View File
@@ -1,8 +1,10 @@
import numpy as np
from skimage.feature import greycomatrix, greycoprops
from skimage.feature._texture import greycomatrix, greycoprops, \
local_binary_pattern, bit_rotate_right
class TestGLCM():
def setup(self):
self.image = np.array([[0, 0, 1, 1],
[0, 0, 1, 1],
@@ -140,5 +142,66 @@ 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=np.uint8)
def test_bit_rotate_right(self):
np.testing.assert_equal(bit_rotate_right(11, 4), 13)
def test_default(self):
lbp = local_binary_pattern(self.image, 8, 1)
ref = np.array([[ 0, 251, 0, 255, 96, 255],
[143, 0, 20, 153, 64, 56],
[238, 255, 12, 191, 0, 252],
[129, 0, 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, 0, 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, 0, 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. , 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()