ENH: Added a gallery example of GLCM textures

This commit is contained in:
Neil Yager
2011-11-21 09:51:20 +00:00
parent cd8524b93d
commit 1661015af1
5 changed files with 183 additions and 66 deletions
+1 -1
View File
@@ -47,7 +47,7 @@ modified to work as part of the scikit, others may be lacking in documentation
or tests.
* :strike:`Connected components`
* `Grey-level co-occurrence matrices <http://mentat.za.net/hg>`_
* :strike:`Grey-level co-occurrence matrices`
* Marching squares
* Nadav's bilateral filtering (first compare against CellProfiler's
code, based on http://groups.csail.mit.edu/graphics/bilagrid/bilagrid_web.pdf)
+87
View File
@@ -0,0 +1,87 @@
"""
=====================
GLCM Texture Features
=====================
This module provides an example of texture classification using grey
level co-occurance matrices (GLCMs). A GLCM is a histogram of
co-occuring greyscale values at a given offset over an image.
In this example, samples of two different textures are extracted from
an image: grassy areas and sky areas. For each patch, a GLCM with
a horizontal offset of 5 is computed. Next, two features of the
GLCM matrices are computed: dissimilarity and correlation. These are
plotted to illustrate that the classes form clusters in feature space.
In a typical classification problem, the final step (not included in
this example) would be to train a classifier, such as logistic
regression, to label image patches from new images.
"""
import os
from skimage.feature import compute_glcm, compute_glcm_prop
from skimage.io import imread
from skimage import data_dir
import matplotlib.pyplot as plt
PATCH_SIZE = 21
# open the camera image
image = imread(os.path.join(data_dir, 'camera.png'))
if False:
plt.figure()
plt.imshow(image)
plt.show()
import sys
sys.exit()
# select some patches from grassy areas of the image
locations = [(474, 291), (440, 433), (466, 18), (462, 236)]
grass_patches = []
for loc in locations:
grass_patches.append(image[loc[0]:loc[0] + PATCH_SIZE,
loc[1]:loc[1] + PATCH_SIZE])
# select some patches from sky areas of the image
locations = [(54, 48), (21, 233), (90, 380), (195, 330)]
sky_patches = []
for loc in locations:
sky_patches.append(image[loc[0]:loc[0] + PATCH_SIZE,
loc[1]:loc[1] + PATCH_SIZE])
# compute some GLCM properties each patch
xs = []
ys = []
for i, patch in enumerate(grass_patches + sky_patches):
glcm = compute_glcm(patch, [5], [0], 256, symmetric=True, normed=True)
xs.append(compute_glcm_prop(glcm, 'dissimilarity')[0, 0])
ys.append(compute_glcm_prop(glcm, 'correlation')[0, 0])
# display the image patches
plt.figure(figsize=(8, 8))
for i, patch in enumerate(grass_patches):
plt.subplot(3, len(grass_patches), i+1)
plt.imshow(patch, cmap=plt.cm.gray, interpolation='nearest',
vmin=0, vmax=255)
plt.xlabel('Grass %d'%(i + 1))
for i, patch in enumerate(sky_patches):
plt.subplot(3, len(grass_patches), i+len(grass_patches)+1)
plt.imshow(patch, cmap=plt.cm.gray, interpolation='nearest',
vmin=0, vmax=255)
plt.xlabel('Sky %d'%(i + 1))
# for each patch, plot (dissimilarity, correlation)
plt.subplot(3, 1, 3)
plt.plot(xs[:len(grass_patches)], ys[:len(grass_patches)], 'go',
label='Grass')
plt.plot(xs[len(grass_patches):], ys[len(grass_patches):], 'bo',
label='Sky')
plt.xlabel('GLCM Dissimilarity')
plt.ylabel('GLVM Correlation')
plt.legend()
# display the patches and plot
plt.suptitle('Grey level co-occurance matrix features', fontsize=14)
plt.show()
+1 -1
View File
@@ -1,2 +1,2 @@
from hog import hog
from greycomatrix import glcm, compute_glcm_prop
from greycomatrix import compute_glcm, compute_glcm_prop
+66 -41
View File
@@ -1,50 +1,47 @@
"""
Compute grey level co-occurrence matrices (GLCM) to characterize
image textures.
Compute grey level co-occurrence matrices (GLCM) and associated
properties to characterize image textures.
"""
import numpy as np
import skimage.util
def glcm(image, distances, angles, levels=256, symmetric=False,
def compute_glcm(image, distances, angles, levels=256, symmetric=False,
normed=False):
"""Calculate the grey-level co-occurrence matrix of a grey-level
image.
"""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. It can be used to
extract features from textured areas of an image.
greyscale values at a given offset over an image.
Parameters
----------
image : (M,N) ndarray
Input image. The input image is converted to the uint8 data
type.
distances : (K,) ndarray
Histogram distance offsets
angles : (L,) ndarray
Histogram angles in radians
levels : int
image : ndarray
Input image, which is converted to the uint8 data type.
distances : array_like
List of histogram distance offsets.
angles : array_like
List of histogram 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).
symmetric : bool
(typically 256 for an 8-bit image). The default is 256.
symmetric : bool, optional
If True, the output matrix P 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.
normed : bool
are accumulated when (i, j) is encountered. The default is False.
normed : bool, optional
If True, normalize the result by dividing by the number of
possible outcomes. The elements of the resulting matrix sum
to 1.
to 1. The default is False.
Returns
-------
P : 4-dimensional 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. I
out : 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`.
References
----------
@@ -63,7 +60,7 @@ def glcm(image, distances, angles, levels=256, symmetric=False,
... [0, 0, 1, 1],
... [0, 2, 2, 2],
... [2, 2, 3, 3]], dtype=np.uint8)
>>> result = glcm(image, [1], [0, np.pi/2], 4)
>>> result = compute_glcm(image, [1], [0, np.pi/2], 4)
>>> result[:, :, 0, 0]
array([[2, 2, 1, 0],
[0, 2, 0, 0],
@@ -124,10 +121,42 @@ def glcm(image, distances, angles, levels=256, symmetric=False,
return out
def compute_glcm_prop(glcm, prop='contrast'):
"""Calculate texture properties of a GLCM.
TODO: rest of docstring, including math
Compute a feature of a grey level co-occurance matrix to serve as
a compact summary of the matrix. The properties are computed as
follows:
- 'contrast': :math:`X`
- 'dissimilarity': :math:`X`
- 'homogeneity': :math:`X`
- 'energy': :math:`X`
- 'correlation': :math:`X`
- 'ASM': :math:`X`
Parameters
----------
glcm : ndarray
Input array. `glcm` 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 : 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
--------
@@ -138,7 +167,8 @@ def compute_glcm_prop(glcm, prop='contrast'):
... [0, 0, 1, 1],
... [0, 2, 2, 2],
... [2, 2, 3, 3]], dtype=np.uint8)
>>> g = glcm(image, [1, 2], [0, np.pi/2], 4, normed=True, symmetric=True)
>>> g = compute_glcm(image, [1, 2], [0, np.pi/2], 4, normed=True,
... symmetric=True)
>>> contrast = compute_glcm_prop(g, 'contrast')
>>> contrast
array([[ 0.58333333, 1. ],
@@ -156,33 +186,33 @@ def compute_glcm_prop(glcm, prop='contrast'):
r = range(num_level)
I, J = np.meshgrid(r, r)
if prop == 'contrast':
weights = (I-J)**2
weights = (I - J) ** 2
elif prop == 'dissimilarity':
weights = np.abs(I-J)
weights = np.abs(I - J)
elif prop == 'homogeneity':
weights = 1./(1.+(I-J)**2)
weights = 1. / (1. + (I - J) ** 2)
elif prop in ['ASM', 'energy', 'correlation']:
pass
else:
raise ValueError('%s is an invalid property'%(prop))
raise ValueError('%s is an invalid property' % (prop))
# compute property for each GLCM
results = np.zeros((num_dist, num_angle), dtype=np.float64)
for d in range(num_dist):
for a in range(num_angle):
if prop == 'energy':
asm = (glcm[:, :, d, a]**2).sum()
asm = (glcm[:, :, d, a] ** 2).sum()
results[d, a] = np.sqrt(asm)
elif prop == 'ASM':
results[d, a] = (glcm[:, :, d, a]**2).sum()
results[d, a] = (glcm[:, :, d, a] ** 2).sum()
elif prop == 'correlation':
g = glcm[:, :, d, a]
mean_i = (I * g).sum()
mean_j = (J * g).sum()
diff_i = I - mean_i
diff_j = J - mean_j
std_i = np.sqrt((g * (diff_i)**2).sum())
std_j = np.sqrt((g * (diff_j)**2).sum())
std_i = np.sqrt((g * (diff_i) ** 2).sum())
std_j = np.sqrt((g * (diff_j) ** 2).sum())
cov = (g * (diff_i * diff_j)).sum()
if std_i < 1e-15 or std_j < 1e-15:
corr = 1.
@@ -194,8 +224,3 @@ def compute_glcm_prop(glcm, prop='contrast'):
results[d, a] = (glcm[:, :, d, a] * weights).sum()
return results
if __name__ == "__main__":
import doctest
doctest.testmod()
+28 -23
View File
@@ -1,5 +1,5 @@
import numpy as np
from skimage.feature import glcm, compute_glcm_prop
from skimage.feature import compute_glcm, compute_glcm_prop
class TestGLCM():
def setup(self):
@@ -9,7 +9,7 @@ class TestGLCM():
[2, 2, 3, 3]], dtype=np.uint8)
def test_output_angles(self):
result = glcm(self.image, [1], [0, np.pi/2], 4)
result = compute_glcm(self.image, [1], [0, np.pi/2], 4)
assert result.shape == (4, 4, 1, 2)
expected1 = np.array([[2, 2, 1, 0],
[0, 2, 0, 0],
@@ -23,7 +23,8 @@ class TestGLCM():
np.testing.assert_array_equal(result[:, :, 0, 1], expected2)
def test_output_symmetric_1(self):
result = glcm(self.image, [1], [np.pi/2], 4, symmetric=True)
result = compute_glcm(self.image, [1], [np.pi/2], 4,
symmetric=True)
assert result.shape == (4, 4, 1, 1)
expected = np.array([[6, 0, 2, 0],
[0, 4, 2, 0],
@@ -32,7 +33,8 @@ class TestGLCM():
np.testing.assert_array_equal(result[:, :, 0, 0], expected)
def test_result_symmetric_2(self):
result = glcm(self.image, [1], [0], 4, symmetric=True)[:, :, 0, 0]
result = compute_glcm(self.image, [1], [0], 4,
symmetric=True)[:, :, 0, 0]
np.testing.assert_array_equal(result, result.transpose())
def test_output_distance(self):
@@ -40,7 +42,7 @@ class TestGLCM():
[1, 0, 0, 1],
[2, 0, 0, 2],
[3, 0, 0, 3]], dtype=np.uint8)
result = glcm(im, [3], [0], 4, symmetric=False)
result = compute_glcm(im, [3], [0], 4, symmetric=False)
expected = np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
@@ -52,7 +54,7 @@ class TestGLCM():
[1],
[2],
[3]], dtype=np.uint8)
result = glcm(im, [1, 2], [0, np.pi/2], 4)
result = compute_glcm(im, [1, 2], [0, np.pi/2], 4)
assert result.shape == (4, 4, 2, 2)
z = np.zeros((4, 4), dtype=np.uint32)
@@ -71,68 +73,71 @@ class TestGLCM():
np.testing.assert_array_equal(result[:, :, 1, 1], e2)
def test_output_empty(self):
result = glcm(self.image, [10], [0], 4)
result = compute_glcm(self.image, [10], [0], 4)
np.testing.assert_array_equal(result[:, :, 0, 0],
np.zeros((4, 4), dtype=np.uint32))
result = glcm(self.image, [10], [0], 4, normed=True)
result = compute_glcm(self.image, [10], [0], 4, normed=True)
np.testing.assert_array_equal(result[:, :, 0, 0],
np.zeros((4, 4), dtype=np.uint32))
def test_normed(self):
result = glcm(self.image, [1, 2, 3], [0, np.pi/2, np.pi], 4,
normed=True)
result = compute_glcm(self.image, [1, 2, 3],
[0, np.pi/2, np.pi], 4, normed=True)
for d in range(result.shape[2]):
for a in range(result.shape[3]):
np.testing.assert_almost_equal(result[:, :, d, a].sum(),
1.0)
def test_contrast(self):
result = glcm(self.image, [1], [0], 4,
normed=True, symmetric=True)
result = compute_glcm(self.image, [1], [0], 4,
normed=True, symmetric=True)
result = np.round(result, 3)
contrast = compute_glcm_prop(result, 'contrast')
np.testing.assert_almost_equal(contrast[0, 0], 0.586)
def test_dissimilarity(self):
result = glcm(self.image, [1], [0], 4,
normed=True, symmetric=True)
result = compute_glcm(self.image, [1], [0], 4,
normed=True, symmetric=True)
result = np.round(result, 3)
dissimilarity = compute_glcm_prop(result, 'dissimilarity')
np.testing.assert_almost_equal(dissimilarity[0, 0], 0.418)
def test_dissimilarity_2(self):
result = glcm(self.image, [1], [np.pi/2], 4,
normed=True, symmetric=True)
result = compute_glcm(self.image, [1], [np.pi/2], 4,
normed=True, symmetric=True)
result = np.round(result, 3)
dissimilarity = compute_glcm_prop(result, 'dissimilarity')[0, 0]
np.testing.assert_almost_equal(dissimilarity, 0.664)
def test_invalid_property(self):
result = glcm(self.image, [1], [0], 4)
result = compute_glcm(self.image, [1], [0], 4)
np.testing.assert_raises(ValueError, compute_glcm_prop,
result, 'ABC')
def test_homogeneity(self):
result = glcm(self.image, [1], [0], 4, normed=True, symmetric=True)
result = compute_glcm(self.image, [1], [0], 4, normed=True,
symmetric=True)
homogeneity = compute_glcm_prop(result, 'homogeneity')[0, 0]
np.testing.assert_almost_equal(homogeneity, 0.80833333)
def test_energy(self):
result = glcm(self.image, [1], [0], 4, normed=True, symmetric=True)
result = compute_glcm(self.image, [1], [0], 4, normed=True,
symmetric=True)
energy = compute_glcm_prop(result, 'energy')[0, 0]
np.testing.assert_almost_equal(energy, 0.38188131)
def test_correlation(self):
result = glcm(self.image, [1], [0], 4, normed=True, symmetric=True)
result = compute_glcm(self.image, [1], [0], 4, normed=True,
symmetric=True)
energy = compute_glcm_prop(result, 'correlation')[0, 0]
np.testing.assert_almost_equal(energy, 0.71953255)
def test_uniform_properties(self):
im = np.ones((4, 4), dtype=np.uint8)
result = glcm(im, [1, 2], [0, np.pi/2], 4, normed=True,
symmetric=True)
result = compute_glcm(im, [1, 2], [0, np.pi/2], 4, normed=True,
symmetric=True)
for prop in ['contrast', 'dissimilarity', 'homogeneity',
'energy', 'correlation']:
'energy', 'correlation', 'ASM']:
compute_glcm_prop(result, prop)
if __name__ == '__main__':