ENH: Updates based on review comments

This commit is contained in:
Neil Yager
2011-11-21 09:51:20 +00:00
parent 1d540c2b56
commit ab780b1ecd
5 changed files with 133 additions and 121 deletions
+1 -1
View File
@@ -79,7 +79,7 @@
Windows packaging and Python 3 compatibility.
- Neil Yager
Skeletonization.
Skeletonization and grey level co-occurrence matrices.
- Nelle Varoquaux
Renaming of the package to ``skimage``.
+46 -39
View File
@@ -3,80 +3,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.
This example illustrates texture classification using texture
classification using grey level co-occurrence matrices (GLCMs).
A GLCM is a histogram of co-occurring 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
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.
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
from skimage.feature import greycomatrix, greycoprops
from skimage import data
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()
image = data.camera()
# select some patches from grassy areas of the image
locations = [(474, 291), (440, 433), (466, 18), (462, 236)]
grass_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,
for loc in grass_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_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,
for loc in sky_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])
glcm = greycomatrix(patch, [5], [0], 256, symmetric=True, normed=True)
xs.append(greycoprops(glcm, 'dissimilarity')[0, 0])
ys.append(greycoprops(glcm, 'correlation')[0, 0])
# create the figure
plt.figure(figsize=(8, 8))
# 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',
plt.subplot(3, len(grass_patches), len(grass_patches) * 1 + i + 1)
plt.imshow(patch, cmap=plt.cm.gray, interpolation='nearest',
vmin=0, vmax=255)
plt.xlabel('Grass %d'%(i + 1))
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))
plt.subplot(3, len(grass_patches), len(grass_patches) * 2 + i + 1)
plt.imshow(patch, cmap=plt.cm.gray, interpolation='nearest',
vmin=0, vmax=255)
plt.xlabel('Sky %d' % (i + 1))
# display original image with locations of patches
plt.subplot(3, 2, 1)
plt.imshow(image, cmap=plt.cm.gray, interpolation='nearest',
vmin=0, vmax=255)
for (y, x) in grass_locations:
plt.plot(x + PATCH_SIZE / 2, y + PATCH_SIZE / 2, 'gs')
for (y, x) in sky_locations:
plt.plot(x + PATCH_SIZE / 2, y + PATCH_SIZE / 2, 'bs')
plt.xlabel('Original Image')
plt.xticks([])
plt.yticks([])
# for each patch, plot (dissimilarity, correlation)
plt.subplot(3, 1, 3)
plt.plot(xs[:len(grass_patches)], ys[:len(grass_patches)], 'go',
plt.subplot(3, 2, 2)
plt.plot(xs[:len(grass_patches)], ys[:len(grass_patches)], 'go',
label='Grass')
plt.plot(xs[len(grass_patches):], ys[len(grass_patches):], 'bo',
plt.plot(xs[len(grass_patches):], ys[len(grass_patches):], 'bo',
label='Sky')
plt.xlabel('GLCM Dissimilarity')
plt.ylabel('GLVM Correlation')
+1 -1
View File
@@ -1,2 +1,2 @@
from hog import hog
from greycomatrix import compute_glcm, compute_glcm_prop
from greycomatrix import greycomatrix, greycoprops
+51 -45
View File
@@ -9,8 +9,8 @@ import skimage.util
from _greycomatrix import _glcm_loop
def compute_glcm(image, distances, angles, levels=256, symmetric=False,
normed=False):
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
@@ -19,7 +19,8 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False,
Parameters
----------
image : ndarray
Input image, which is converted to the uint8 data type.
Input image. The image is converted to the uint8 data type, so
its range of the image is [0, 255].
distances : array_like
List of pixel pair distance offsets.
angles : array_like
@@ -27,19 +28,21 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False,
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 default is 256.
(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. The default is False.
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 the result by dividing by the number of
possible outcomes. The elements of the resulting matrix sum
to 1. The default is False.
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
-------
hist : ndarray
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
@@ -52,18 +55,19 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False,
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 = compute_glcm(image, [1], [0, np.pi/2], 4)
>>> result = greycomatrix(image, [1], [0, np.pi/2], levels=4)
>>> result[:, :, 0, 0]
array([[2, 2, 1, 0],
[0, 2, 0, 0],
@@ -85,33 +89,29 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False,
assert distances.ndim == 1
assert angles.ndim == 1
hist = np.zeros((levels, levels, len(distances), len(angles)),
dtype=np.uint32, order='C')
P = np.zeros((levels, levels, len(distances), len(angles)),
dtype=np.uint32, order='C')
# count co-occurances
_glcm_loop(image, distances, angles, levels, hist)
# count co-occurences
_glcm_loop(image, distances, angles, levels, P)
# make each GLMC symmetric
if symmetric:
for d in range(len(distances)):
for a in range(len(angles)):
hist[:, :, d, a] += hist[:, :, d, a].transpose()
P += np.transpose(P, (1, 0, 2, 3))
# normalize each GLMC
if normed:
hist = hist.astype(np.float64)
for d in range(len(distances)):
for a in range(len(angles)):
if np.any(hist[:, :, d, a]):
hist[:, :, d, a] /= hist[:, :, d, a].sum()
P = P.astype(np.float64)
P /= np.apply_over_axes(np.sum, P, axes=(0, 1))
P = np.nan_to_num(P)
return hist
return P
def compute_glcm_prop(P, prop='contrast'):
def greycoprops(P, prop='contrast'):
"""Calculate texture properties of a GLCM.
Compute a feature of a grey level co-occurance matrix to serve as
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:
@@ -136,7 +136,7 @@ def compute_glcm_prop(P, prop='contrast'):
Returns
-------
results : ndarray
results : 2-D ndarray
2-dimensional array. `results[d, a]` is the property 'prop' for
the d'th distance and the a'th angle.
@@ -154,8 +154,8 @@ def compute_glcm_prop(P, prop='contrast'):
... [0, 0, 1, 1],
... [0, 2, 2, 2],
... [2, 2, 3, 3]], dtype=np.uint8)
>>> g = compute_glcm(image, [1, 2], [0, np.pi/2], 4, normed=True,
... symmetric=True)
>>> g = greycomatrix(image, [1, 2], [0, np.pi/2], levels=4,
... normed=True, symmetric=True)
>>> contrast = compute_glcm_prop(g, 'contrast')
>>> contrast
array([[ 0.58333333, 1. ],
@@ -170,8 +170,7 @@ def compute_glcm_prop(P, prop='contrast'):
assert num_angle > 0
# create weights for specified property
r = range(num_level)
I, J = np.meshgrid(r, r)
I, J = np.ogrid[0:num_level, 0:num_level]
if prop == 'contrast':
weights = (I - J) ** 2
elif prop == 'dissimilarity':
@@ -182,17 +181,17 @@ def compute_glcm_prop(P, prop='contrast'):
pass
else:
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 = (P[:, :, d, a] ** 2).sum()
results[d, a] = np.sqrt(asm)
elif prop == 'ASM':
results[d, a] = (P[:, :, d, a] ** 2).sum()
elif prop == 'correlation':
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)
for d in range(num_dist):
for a in range(num_angle):
g = P[:, :, d, a]
mean_i = (I * g).sum()
mean_j = (J * g).sum()
@@ -207,7 +206,14 @@ def compute_glcm_prop(P, prop='contrast'):
corr = cov / (std_i * std_j)
results[d, a] = corr
else:
results[d, a] = (P[:, :, d, a] * weights).sum()
results[d, a] = corr
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
if __name__ == "__main__":
import doctest
doctest.testmod()
+34 -35
View File
@@ -1,5 +1,6 @@
import numpy as np
from skimage.feature import compute_glcm, compute_glcm_prop
from skimage.feature import greycomatrix, greycoprops
class TestGLCM():
def setup(self):
@@ -9,7 +10,7 @@ class TestGLCM():
[2, 2, 3, 3]], dtype=np.uint8)
def test_output_angles(self):
result = compute_glcm(self.image, [1], [0, np.pi/2], 4)
result = greycomatrix(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 +24,7 @@ class TestGLCM():
np.testing.assert_array_equal(result[:, :, 0, 1], expected2)
def test_output_symmetric_1(self):
result = compute_glcm(self.image, [1], [np.pi/2], 4,
result = greycomatrix(self.image, [1], [np.pi / 2], 4,
symmetric=True)
assert result.shape == (4, 4, 1, 1)
expected = np.array([[6, 0, 2, 0],
@@ -32,17 +33,12 @@ class TestGLCM():
[0, 0, 2, 0]], dtype=np.uint32)
np.testing.assert_array_equal(result[:, :, 0, 0], expected)
def test_result_symmetric_2(self):
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):
im = np.array([[0, 0, 0, 0],
[1, 0, 0, 1],
[2, 0, 0, 2],
[3, 0, 0, 3]], dtype=np.uint8)
result = compute_glcm(im, [3], [0], 4, symmetric=False)
result = greycomatrix(im, [3], [0], 4, symmetric=False)
expected = np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
@@ -54,7 +50,7 @@ class TestGLCM():
[1],
[2],
[3]], dtype=np.uint8)
result = compute_glcm(im, [1, 2], [0, np.pi/2], 4)
result = greycomatrix(im, [1, 2], [0, np.pi / 2], 4)
assert result.shape == (4, 4, 2, 2)
z = np.zeros((4, 4), dtype=np.uint32)
@@ -73,73 +69,76 @@ class TestGLCM():
np.testing.assert_array_equal(result[:, :, 1, 1], e2)
def test_output_empty(self):
result = compute_glcm(self.image, [10], [0], 4)
result = greycomatrix(self.image, [10], [0], 4)
np.testing.assert_array_equal(result[:, :, 0, 0],
np.zeros((4, 4), dtype=np.uint32))
result = compute_glcm(self.image, [10], [0], 4, normed=True)
result = greycomatrix(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 = compute_glcm(self.image, [1, 2, 3],
[0, np.pi/2, np.pi], 4, normed=True)
def test_normed_symmetric(self):
result = greycomatrix(self.image, [1, 2, 3],
[0, np.pi / 2, np.pi], 4,
normed=True, symmetric=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)
np.testing.assert_array_equal(result[:, :, d, a],
result[:, :, d, a].transpose())
def test_contrast(self):
result = compute_glcm(self.image, [1], [0], 4,
result = greycomatrix(self.image, [1, 2], [0], 4,
normed=True, symmetric=True)
result = np.round(result, 3)
contrast = compute_glcm_prop(result, 'contrast')
contrast = greycoprops(result, 'contrast')
np.testing.assert_almost_equal(contrast[0, 0], 0.586)
def test_dissimilarity(self):
result = compute_glcm(self.image, [1], [0], 4,
result = greycomatrix(self.image, [1], [0, np.pi / 2], 4,
normed=True, symmetric=True)
result = np.round(result, 3)
dissimilarity = compute_glcm_prop(result, 'dissimilarity')
dissimilarity = greycoprops(result, 'dissimilarity')
np.testing.assert_almost_equal(dissimilarity[0, 0], 0.418)
def test_dissimilarity_2(self):
result = compute_glcm(self.image, [1], [np.pi/2], 4,
result = greycomatrix(self.image, [1, 3], [np.pi/2], 4,
normed=True, symmetric=True)
result = np.round(result, 3)
dissimilarity = compute_glcm_prop(result, 'dissimilarity')[0, 0]
dissimilarity = greycoprops(result, 'dissimilarity')[0, 0]
np.testing.assert_almost_equal(dissimilarity, 0.664)
def test_invalid_property(self):
result = compute_glcm(self.image, [1], [0], 4)
np.testing.assert_raises(ValueError, compute_glcm_prop,
result = greycomatrix(self.image, [1], [0], 4)
np.testing.assert_raises(ValueError, greycoprops,
result, 'ABC')
def test_homogeneity(self):
result = compute_glcm(self.image, [1], [0], 4, normed=True,
result = greycomatrix(self.image, [1], [0, 6], 4, normed=True,
symmetric=True)
homogeneity = compute_glcm_prop(result, 'homogeneity')[0, 0]
homogeneity = greycoprops(result, 'homogeneity')[0, 0]
np.testing.assert_almost_equal(homogeneity, 0.80833333)
def test_energy(self):
result = compute_glcm(self.image, [1], [0], 4, normed=True,
result = greycomatrix(self.image, [1], [0, 4], 4, normed=True,
symmetric=True)
energy = compute_glcm_prop(result, 'energy')[0, 0]
energy = greycoprops(result, 'energy')[0, 0]
np.testing.assert_almost_equal(energy, 0.38188131)
def test_correlation(self):
result = compute_glcm(self.image, [1], [0], 4, normed=True,
result = greycomatrix(self.image, [1, 2], [0], 4, normed=True,
symmetric=True)
energy = compute_glcm_prop(result, 'correlation')[0, 0]
np.testing.assert_almost_equal(energy, 0.71953255)
energy = greycoprops(result, 'correlation')
np.testing.assert_almost_equal(energy[0, 0], 0.71953255)
np.testing.assert_almost_equal(energy[1, 0], 0.41176470)
def test_uniform_properties(self):
im = np.ones((4, 4), dtype=np.uint8)
result = compute_glcm(im, [1, 2], [0, np.pi/2], 4, normed=True,
result = greycomatrix(im, [1, 2], [0, np.pi / 2], 4, normed=True,
symmetric=True)
for prop in ['contrast', 'dissimilarity', 'homogeneity',
'energy', 'correlation', 'ASM']:
compute_glcm_prop(result, prop)
greycoprops(result, prop)
if __name__ == '__main__':
np.testing.run_module_suite()