diff --git a/TASKS.txt b/TASKS.txt index a17ecf35..ad5e2f25 100644 --- a/TASKS.txt +++ b/TASKS.txt @@ -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 `_ + * :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) diff --git a/doc/examples/plot_glcm.py b/doc/examples/plot_glcm.py new file mode 100644 index 00000000..1de6a27b --- /dev/null +++ b/doc/examples/plot_glcm.py @@ -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() diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 58e1e072..c3332ecf 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1,2 +1,2 @@ from hog import hog -from greycomatrix import glcm, compute_glcm_prop +from greycomatrix import compute_glcm, compute_glcm_prop diff --git a/skimage/feature/greycomatrix.py b/skimage/feature/greycomatrix.py index 0557fe2f..c4dd2a98 100644 --- a/skimage/feature/greycomatrix.py +++ b/skimage/feature/greycomatrix.py @@ -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() - \ No newline at end of file diff --git a/skimage/feature/tests/test_glcm.py b/skimage/feature/tests/test_glcm.py index 16c9bb0a..55950d35 100644 --- a/skimage/feature/tests/test_glcm.py +++ b/skimage/feature/tests/test_glcm.py @@ -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__':