From ab780b1ecd9a59a168aab320477e297a50acf77c Mon Sep 17 00:00:00 2001 From: Neil Yager Date: Mon, 14 Nov 2011 14:32:54 +0000 Subject: [PATCH] ENH: Updates based on review comments --- CONTRIBUTORS.txt | 2 +- doc/examples/plot_glcm.py | 85 ++++++++++++++------------ skimage/feature/__init__.py | 2 +- skimage/feature/greycomatrix.py | 96 ++++++++++++++++-------------- skimage/feature/tests/test_glcm.py | 69 +++++++++++---------- 5 files changed, 133 insertions(+), 121 deletions(-) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index bff86ccc..c0816717 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -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``. diff --git a/doc/examples/plot_glcm.py b/doc/examples/plot_glcm.py index 1de6a27b..e30f10d7 100644 --- a/doc/examples/plot_glcm.py +++ b/doc/examples/plot_glcm.py @@ -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') diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index c3332ecf..6b3b7014 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1,2 +1,2 @@ from hog import hog -from greycomatrix import compute_glcm, compute_glcm_prop +from greycomatrix import greycomatrix, greycoprops diff --git a/skimage/feature/greycomatrix.py b/skimage/feature/greycomatrix.py index 595c06fb..00929ee8 100644 --- a/skimage/feature/greycomatrix.py +++ b/skimage/feature/greycomatrix.py @@ -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() diff --git a/skimage/feature/tests/test_glcm.py b/skimage/feature/tests/test_glcm.py index 55950d35..26c4a1f0 100644 --- a/skimage/feature/tests/test_glcm.py +++ b/skimage/feature/tests/test_glcm.py @@ -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() - \ No newline at end of file