From 4fe7965f694f65a982a693ba86bfa4b104eee660 Mon Sep 17 00:00:00 2001 From: Neil Yager Date: Tue, 1 Nov 2011 15:22:00 +0000 Subject: [PATCH 01/13] NF: Initial commit of greylevel co-occurance (pure python) --- skimage/feature/__init__.py | 3 +- skimage/feature/greycomatrix.py | 112 +++++++++++++++++++++++++++++ skimage/feature/tests/test_glcm.py | 90 +++++++++++++++++++++++ 3 files changed, 204 insertions(+), 1 deletion(-) create mode 100644 skimage/feature/greycomatrix.py create mode 100644 skimage/feature/tests/test_glcm.py diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index eb6faa62..4a0b2955 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1 +1,2 @@ -from hog import hog \ No newline at end of file +from hog import hog +from greycomatrix import glcm diff --git a/skimage/feature/greycomatrix.py b/skimage/feature/greycomatrix.py new file mode 100644 index 00000000..82d8f86b --- /dev/null +++ b/skimage/feature/greycomatrix.py @@ -0,0 +1,112 @@ +""" +Compute grey level co-occurrence matrices (GLCM) to characterize +image textures. +""" + +import numpy as np +import skimage.util + + +def glcm(image, distances, angles, levels=256, symmetric=False, + normal=False): + """Calculate the grey-level co-occurrence matrix of a grey-level + image. + + 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. + + 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 + 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 + 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. + normal : bool + If True, normalize the result by dividing by the number of + possible outcomes + + 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. + + 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 = glcm(image, [1], [0, np.pi/2], 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) + + """ + image = skimage.util.img_as_ubyte(image) + assert image.ndim == 2 + assert image.min() >= 0 + assert image.max() < levels + distances = np.asarray(distances) + angles = np.asarray(angles) + assert distances.ndim == 1 + assert angles.ndim == 1 + + rows, cols = image.shape + out = np.zeros((levels, levels, len(distances), len(angles)), + dtype=np.uint32) + + for a_idx, angle in enumerate(angles): + for d_idx, distance in enumerate(distances): + for r in range(rows): + for c in range(cols): + i = image[r, c] + + # compute the location of the offset pixel + row = r + int(np.round(np.sin(angle) * distance)) + col = c + int(np.round(np.cos(angle) * distance)) + + # make sure the offset is within bounds + if row >= 0 and row < rows and \ + col >= 0 and col < cols: + j = image[row, col] + + if i >= 0 and i < levels and \ + j >= 0 and j < levels: + out[i, j, d_idx, a_idx] += 1 + if symmetric: + out[j, i, d_idx, a_idx] += 1 + + # normalize + if normal: + out = out.astype(np.float64) / out.sum() + + return out + +if __name__ == "__main__": + import doctest + doctest.testmod() diff --git a/skimage/feature/tests/test_glcm.py b/skimage/feature/tests/test_glcm.py new file mode 100644 index 00000000..58a623c4 --- /dev/null +++ b/skimage/feature/tests/test_glcm.py @@ -0,0 +1,90 @@ +import numpy as np +from skimage.feature import glcm + +class TestGLCM(): + def test_output_angles(self): + image = np.array([[0, 0, 1, 1], + [0, 0, 1, 1], + [0, 2, 2, 2], + [2, 2, 3, 3]], dtype=np.uint8) + result = glcm(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], + [0, 0, 3, 1], + [0, 0, 0, 1]], dtype=np.uint32) + np.testing.assert_array_equal(result[:, :, 0, 0], expected1) + expected2 = np.array([[3, 0, 2, 0], + [0, 2, 2, 0], + [0, 0, 1, 2], + [0, 0, 0, 0]], dtype=np.uint32) + np.testing.assert_array_equal(result[:, :, 0, 1], expected2) + + def test_output_symmetric_1(self): + image = np.array([[0, 0, 1, 1], + [0, 0, 1, 1], + [0, 2, 2, 2], + [2, 2, 3, 3]], dtype=np.uint8) + result = glcm(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], + [2, 2, 2, 2], + [0, 0, 2, 0]], dtype=np.uint32) + np.testing.assert_array_equal(result[:, :, 0, 0], expected) + + def test_result_symmetric_2(self): + image = np.array([[0, 0, 1, 1], + [0, 0, 1, 1], + [0, 2, 2, 2], + [2, 2, 3, 3]], dtype=np.uint8) + result = glcm(image, [1], [0], 4, symmetric=True)[:, :, 0, 0] + np.testing.assert_array_equal(result, result.transpose()) + + def test_output_distance(self): + image = np.array([[0, 0, 0, 0], + [1, 0, 0, 1], + [2, 0, 0, 2], + [3, 0, 0, 3]], dtype=np.uint8) + result = glcm(image, [3], [0], 4, symmetric=False) + expected = np.array([[1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 1, 0], + [0, 0, 0, 1]], dtype=np.uint32) + np.testing.assert_array_equal(result[:, :, 0, 0], expected) + + def test_output_combo(self): + image = np.array([[0], + [1], + [2], + [3]], dtype=np.uint8) + result = glcm(image, [1, 2], [0, np.pi/2], 4) + assert result.shape == (4, 4, 2, 2) + + z = np.zeros((4, 4), dtype=np.uint32) + e1 = np.array([[0, 1, 0, 0], + [0, 0, 1, 0], + [0, 0, 0, 1], + [0, 0, 0, 0]], dtype=np.uint32) + e2 = np.array([[0, 0, 1, 0], + [0, 0, 0, 1], + [0, 0, 0, 0], + [0, 0, 0, 0]], dtype=np.uint32) + + np.testing.assert_array_equal(result[:, :, 0, 0], z) + np.testing.assert_array_equal(result[:, :, 1, 0], z) + np.testing.assert_array_equal(result[:, :, 0, 1], e1) + np.testing.assert_array_equal(result[:, :, 1, 1], e2) + + def test_normed(self): + image = np.array([[0, 0, 1, 1], + [0, 0, 1, 1], + [0, 2, 2, 2], + [2, 2, 3, 3]], dtype=np.uint8) + result = glcm(image, [1], [0], 4, normal=True) + np.testing.assert_almost_equal(result.sum(), 1.0) + + +if __name__ == '__main__': + np.testing.run_module_suite() + \ No newline at end of file From cd8524b93db803336f42cb7853fdbab0345a1a19 Mon Sep 17 00:00:00 2001 From: Neil Yager Date: Wed, 2 Nov 2011 12:23:27 +0000 Subject: [PATCH 02/13] ENH: Compute some standard propertices for GLCMs --- skimage/feature/__init__.py | 2 +- skimage/feature/greycomatrix.py | 107 ++++++++++++++++++++++--- skimage/feature/tests/test_glcm.py | 122 ++++++++++++++++++++--------- 3 files changed, 185 insertions(+), 46 deletions(-) diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 4a0b2955..58e1e072 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -1,2 +1,2 @@ from hog import hog -from greycomatrix import glcm +from greycomatrix import glcm, compute_glcm_prop diff --git a/skimage/feature/greycomatrix.py b/skimage/feature/greycomatrix.py index 82d8f86b..0557fe2f 100644 --- a/skimage/feature/greycomatrix.py +++ b/skimage/feature/greycomatrix.py @@ -8,7 +8,7 @@ import skimage.util def glcm(image, distances, angles, levels=256, symmetric=False, - normal=False): + normed=False): """Calculate the grey-level co-occurrence matrix of a grey-level image. @@ -33,9 +33,10 @@ def glcm(image, distances, angles, levels=256, symmetric=False, 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. - normal : bool + normed : bool If True, normalize the result by dividing by the number of - possible outcomes + possible outcomes. The elements of the resulting matrix sum + to 1. Returns ------- @@ -43,7 +44,15 @@ def glcm(image, distances, angles, levels=256, symmetric=False, 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. + grey-level i. I + + References + ---------- + .. [1] The GLCM Tutorial Home Page, + http://www.fp.ucalgary.ca/mhallbey/tutorial.htm + .. [2] Pattern Recognition Engineering, Morton Nadler & Eric P. + Smith + Examples -------- @@ -98,15 +107,95 @@ def glcm(image, distances, angles, levels=256, symmetric=False, if i >= 0 and i < levels and \ j >= 0 and j < levels: out[i, j, d_idx, a_idx] += 1 - if symmetric: - out[j, i, d_idx, a_idx] += 1 - # normalize - if normal: - out = out.astype(np.float64) / out.sum() + # make each GLMC symmetric + if symmetric: + for d in range(len(distances)): + for a in range(len(angles)): + out[:, :, d, a] += out[:, :, d, a].transpose() + + # normalize each GLMC individually + if normed: + out = out.astype(np.float64) + for d in range(len(distances)): + for a in range(len(angles)): + if np.any(out[:, :, d, a]): + out[:, :, d, a] /= out[:, :, d, a].sum() return out +def compute_glcm_prop(glcm, prop='contrast'): + """Calculate texture properties of a GLCM. + + TODO: rest of docstring, including math + + 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 = glcm(image, [1, 2], [0, np.pi/2], 4, normed=True, symmetric=True) + >>> contrast = compute_glcm_prop(g, 'contrast') + >>> contrast + array([[ 0.58333333, 1. ], + [ 1.25 , 2.75 ]]) + + """ + + assert glcm.ndim == 4 + (num_level, num_level2, num_dist, num_angle) = glcm.shape + assert num_level == num_level2 + assert num_dist > 0 + assert num_angle > 0 + + # create weights for specified property + r = range(num_level) + I, J = np.meshgrid(r, r) + 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 + 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() + results[d, a] = np.sqrt(asm) + elif prop == 'ASM': + 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()) + cov = (g * (diff_i * diff_j)).sum() + if std_i < 1e-15 or std_j < 1e-15: + corr = 1. + else: + corr = cov / (std_i * std_j) + + results[d, a] = corr + else: + 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 58a623c4..16c9bb0a 100644 --- a/skimage/feature/tests/test_glcm.py +++ b/skimage/feature/tests/test_glcm.py @@ -1,13 +1,15 @@ import numpy as np -from skimage.feature import glcm +from skimage.feature import glcm, compute_glcm_prop class TestGLCM(): + def setup(self): + self.image = np.array([[0, 0, 1, 1], + [0, 0, 1, 1], + [0, 2, 2, 2], + [2, 2, 3, 3]], dtype=np.uint8) + def test_output_angles(self): - image = np.array([[0, 0, 1, 1], - [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 = 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], @@ -20,12 +22,8 @@ class TestGLCM(): [0, 0, 0, 0]], dtype=np.uint32) np.testing.assert_array_equal(result[:, :, 0, 1], expected2) - def test_output_symmetric_1(self): - image = np.array([[0, 0, 1, 1], - [0, 0, 1, 1], - [0, 2, 2, 2], - [2, 2, 3, 3]], dtype=np.uint8) - result = glcm(image, [1], [np.pi/2], 4, symmetric=True) + def test_output_symmetric_1(self): + result = 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], @@ -33,32 +31,28 @@ class TestGLCM(): [0, 0, 2, 0]], dtype=np.uint32) np.testing.assert_array_equal(result[:, :, 0, 0], expected) - def test_result_symmetric_2(self): - image = np.array([[0, 0, 1, 1], - [0, 0, 1, 1], - [0, 2, 2, 2], - [2, 2, 3, 3]], dtype=np.uint8) - result = glcm(image, [1], [0], 4, symmetric=True)[:, :, 0, 0] + def test_result_symmetric_2(self): + result = glcm(self.image, [1], [0], 4, symmetric=True)[:, :, 0, 0] np.testing.assert_array_equal(result, result.transpose()) def test_output_distance(self): - image = np.array([[0, 0, 0, 0], - [1, 0, 0, 1], - [2, 0, 0, 2], - [3, 0, 0, 3]], dtype=np.uint8) - result = glcm(image, [3], [0], 4, symmetric=False) + im = np.array([[0, 0, 0, 0], + [1, 0, 0, 1], + [2, 0, 0, 2], + [3, 0, 0, 3]], dtype=np.uint8) + result = glcm(im, [3], [0], 4, symmetric=False) expected = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]], dtype=np.uint32) - np.testing.assert_array_equal(result[:, :, 0, 0], expected) + np.testing.assert_array_equal(result[:, :, 0, 0], expected) def test_output_combo(self): - image = np.array([[0], - [1], - [2], - [3]], dtype=np.uint8) - result = glcm(image, [1, 2], [0, np.pi/2], 4) + im = np.array([[0], + [1], + [2], + [3]], dtype=np.uint8) + result = glcm(im, [1, 2], [0, np.pi/2], 4) assert result.shape == (4, 4, 2, 2) z = np.zeros((4, 4), dtype=np.uint32) @@ -76,15 +70,71 @@ class TestGLCM(): np.testing.assert_array_equal(result[:, :, 0, 1], e1) np.testing.assert_array_equal(result[:, :, 1, 1], e2) - def test_normed(self): - image = np.array([[0, 0, 1, 1], - [0, 0, 1, 1], - [0, 2, 2, 2], - [2, 2, 3, 3]], dtype=np.uint8) - result = glcm(image, [1], [0], 4, normal=True) - np.testing.assert_almost_equal(result.sum(), 1.0) + def test_output_empty(self): + result = 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) + 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) + 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 = 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 = 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 = 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) + 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) + 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) + 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) + 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) + for prop in ['contrast', 'dissimilarity', 'homogeneity', + 'energy', 'correlation']: + compute_glcm_prop(result, prop) + if __name__ == '__main__': np.testing.run_module_suite() \ No newline at end of file From 1661015af17c59412de655aa2157657394b01d77 Mon Sep 17 00:00:00 2001 From: Neil Yager Date: Wed, 2 Nov 2011 16:09:31 +0000 Subject: [PATCH 03/13] ENH: Added a gallery example of GLCM textures --- TASKS.txt | 2 +- doc/examples/plot_glcm.py | 87 +++++++++++++++++++++++ skimage/feature/__init__.py | 2 +- skimage/feature/greycomatrix.py | 107 ++++++++++++++++++----------- skimage/feature/tests/test_glcm.py | 51 +++++++------- 5 files changed, 183 insertions(+), 66 deletions(-) create mode 100644 doc/examples/plot_glcm.py 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__': From 59ff9218b4ec342ace8a7c7c81a0896edbc2d3cd Mon Sep 17 00:00:00 2001 From: Neil Yager Date: Thu, 3 Nov 2011 08:41:08 +0000 Subject: [PATCH 04/13] DOC: Added latex formulas to docstring --- skimage/feature/greycomatrix.py | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/skimage/feature/greycomatrix.py b/skimage/feature/greycomatrix.py index c4dd2a98..98631e28 100644 --- a/skimage/feature/greycomatrix.py +++ b/skimage/feature/greycomatrix.py @@ -122,24 +122,25 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False, return out -def compute_glcm_prop(glcm, prop='contrast'): +def compute_glcm_prop(P, prop='contrast'): """Calculate texture properties of a GLCM. 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` + + - '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}\\left|i-j\\right|` + - '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 ---------- - glcm : ndarray - Input array. `glcm` is the grey-level co-occurrence histogram + 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 @@ -176,8 +177,8 @@ def compute_glcm_prop(glcm, prop='contrast'): """ - assert glcm.ndim == 4 - (num_level, num_level2, num_dist, num_angle) = glcm.shape + 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 @@ -201,12 +202,12 @@ def compute_glcm_prop(glcm, prop='contrast'): for d in range(num_dist): for a in range(num_angle): if prop == 'energy': - asm = (glcm[:, :, d, a] ** 2).sum() + asm = (P[:, :, d, a] ** 2).sum() results[d, a] = np.sqrt(asm) elif prop == 'ASM': - results[d, a] = (glcm[:, :, d, a] ** 2).sum() + results[d, a] = (P[:, :, d, a] ** 2).sum() elif prop == 'correlation': - g = glcm[:, :, d, a] + g = P[:, :, d, a] mean_i = (I * g).sum() mean_j = (J * g).sum() diff_i = I - mean_i @@ -221,6 +222,6 @@ def compute_glcm_prop(glcm, prop='contrast'): results[d, a] = corr else: - results[d, a] = (glcm[:, :, d, a] * weights).sum() + results[d, a] = (P[:, :, d, a] * weights).sum() return results From 80a83c182f75be17c57c27c7628123a27cd47cb2 Mon Sep 17 00:00:00 2001 From: Neil Yager Date: Thu, 3 Nov 2011 17:02:42 +0000 Subject: [PATCH 05/13] ENH: Use Cython to compute GLCM --- skimage/feature/_greycomatrix.pyx | 67 +++++++++++++++++++++++++++++++ skimage/feature/greycomatrix.py | 32 +++++---------- skimage/feature/setup.py | 30 ++++++++++++++ 3 files changed, 106 insertions(+), 23 deletions(-) create mode 100644 skimage/feature/_greycomatrix.pyx create mode 100644 skimage/feature/setup.py diff --git a/skimage/feature/_greycomatrix.pyx b/skimage/feature/_greycomatrix.pyx new file mode 100644 index 00000000..b9ff7f23 --- /dev/null +++ b/skimage/feature/_greycomatrix.pyx @@ -0,0 +1,67 @@ +"""Cython implementation for computing a grey level co-occurance matrix +""" + +import numpy as np +cimport numpy as np +cimport cython + +cdef extern from "math.h": + double sin(double) + double cos(double) + +@cython.boundscheck(False) +def _glcm_loop(np.ndarray[dtype=np.uint8_t, ndim=2, + negative_indices=False, mode='c'] image, + np.ndarray[dtype=np.float64_t, ndim=1, + negative_indices=False, mode='c'] distances, + np.ndarray[dtype=np.float64_t, ndim=1, + negative_indices=False, mode='c'] angles, + int levels, + np.ndarray[dtype=np.uint32_t, ndim=4, + negative_indices=False, mode='c'] out + ): + """Perform co-occurnace matrix accumulation + + Parameters + ---------- + image : ndarray + Input image, which is converted to the uint8 data type. + distances : ndarray + List of pixel pair distance offsets. + angles : ndarray + List of pixel pair angles in radians. + levels : int + 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) + out : ndarray + On input a 4D array of zeros, and on output it contains + the results of the GLCM computation. + + """ + cdef: + np.int32_t a_inx, d_idx + np.int32_t r, c, rows, cols, row, col + np.int32_t i, j + + rows = image.shape[0] + cols = image.shape[1] + + for a_idx, angle in enumerate(angles): + for d_idx, distance in enumerate(distances): + for r in range(rows): + for c in range(cols): + i = image[r, c] + + # compute the location of the offset pixel + row = r + (sin(angle) * distance + 0.5) + col = c + (cos(angle) * distance + 0.5); + + # make sure the offset is within bounds + if row >= 0 and row < rows and \ + col >= 0 and col < cols: + j = image[row, col] + + if i >= 0 and i < levels and \ + j >= 0 and j < levels: + out[i, j, d_idx, a_idx] += 1 diff --git a/skimage/feature/greycomatrix.py b/skimage/feature/greycomatrix.py index 98631e28..b282dcff 100644 --- a/skimage/feature/greycomatrix.py +++ b/skimage/feature/greycomatrix.py @@ -6,6 +6,8 @@ properties to characterize image textures. import numpy as np import skimage.util +from _greycomatrix import _glcm_loop + def compute_glcm(image, distances, angles, levels=256, symmetric=False, normed=False): @@ -19,9 +21,9 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False, image : ndarray Input image, which is converted to the uint8 data type. distances : array_like - List of histogram distance offsets. + List of pixel pair distance offsets. angles : array_like - List of histogram angles in radians. + 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 @@ -77,8 +79,8 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False, assert image.ndim == 2 assert image.min() >= 0 assert image.max() < levels - distances = np.asarray(distances) - angles = np.asarray(angles) + distances = np.ascontiguousarray(distances, dtype=np.float64) + angles = np.ascontiguousarray(angles, dtype=np.float64) assert distances.ndim == 1 assert angles.ndim == 1 @@ -86,24 +88,8 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False, out = np.zeros((levels, levels, len(distances), len(angles)), dtype=np.uint32) - for a_idx, angle in enumerate(angles): - for d_idx, distance in enumerate(distances): - for r in range(rows): - for c in range(cols): - i = image[r, c] - - # compute the location of the offset pixel - row = r + int(np.round(np.sin(angle) * distance)) - col = c + int(np.round(np.cos(angle) * distance)) - - # make sure the offset is within bounds - if row >= 0 and row < rows and \ - col >= 0 and col < cols: - j = image[row, col] - - if i >= 0 and i < levels and \ - j >= 0 and j < levels: - out[i, j, d_idx, a_idx] += 1 + # count co-occurances + _glcm_loop(image, distances, angles, levels, out) # make each GLMC symmetric if symmetric: @@ -111,7 +97,7 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False, for a in range(len(angles)): out[:, :, d, a] += out[:, :, d, a].transpose() - # normalize each GLMC individually + # normalize each GLMC if normed: out = out.astype(np.float64) for d in range(len(distances)): diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py new file mode 100644 index 00000000..d2ab7d22 --- /dev/null +++ b/skimage/feature/setup.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python + +import os +from skimage._build import cython + +base_path = os.path.abspath(os.path.dirname(__file__)) + +def configuration(parent_package='', top_path=None): + from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs + + config = Configuration('feature', parent_package, top_path) + config.add_data_dir('tests') + + cython(['_greycomatrix.pyx'], working_path=base_path) + + config.add_extension('_greycomatrix', sources=['_greycomatrix.c'], + include_dirs=[get_numpy_include_dirs()]) + + return config + +if __name__ == '__main__': + from numpy.distutils.core import setup + setup(maintainer = 'Scikits-image Developers', + author = 'Scikits-image Developers', + maintainer_email = 'scikits-image@googlegroups.com', + description = 'Features', + url = 'https://github.com/scikits-image/scikits-image', + license = 'SciPy License (BSD Style)', + **(configuration(top_path='').todict()) + ) From 1d540c2b56a0ee44a60734f7f04b8675de82e938 Mon Sep 17 00:00:00 2001 From: Neil Yager Date: Fri, 11 Nov 2011 12:27:48 +0000 Subject: [PATCH 06/13] BUG: Make sure image array is contiguous --- skimage/feature/greycomatrix.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/skimage/feature/greycomatrix.py b/skimage/feature/greycomatrix.py index b282dcff..595c06fb 100644 --- a/skimage/feature/greycomatrix.py +++ b/skimage/feature/greycomatrix.py @@ -1,5 +1,5 @@ """ -Compute grey level co-occurrence matrices (GLCM) and associated +Compute grey level co-occurrence matrices (GLCMs) and associated properties to characterize image textures. """ @@ -39,11 +39,12 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False, Returns ------- - out : ndarray + hist : 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`. + grey-level `i`. If `normed` is `False`, the output is of + type uint32, otherwise it is float64. References ---------- @@ -75,7 +76,7 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False, [0, 0, 0, 0]], dtype=uint32) """ - image = skimage.util.img_as_ubyte(image) + image = np.ascontiguousarray(skimage.util.img_as_ubyte(image)) assert image.ndim == 2 assert image.min() >= 0 assert image.max() < levels @@ -84,28 +85,27 @@ def compute_glcm(image, distances, angles, levels=256, symmetric=False, assert distances.ndim == 1 assert angles.ndim == 1 - rows, cols = image.shape - out = np.zeros((levels, levels, len(distances), len(angles)), - dtype=np.uint32) + hist = np.zeros((levels, levels, len(distances), len(angles)), + dtype=np.uint32, order='C') # count co-occurances - _glcm_loop(image, distances, angles, levels, out) + _glcm_loop(image, distances, angles, levels, hist) # make each GLMC symmetric if symmetric: for d in range(len(distances)): for a in range(len(angles)): - out[:, :, d, a] += out[:, :, d, a].transpose() + hist[:, :, d, a] += hist[:, :, d, a].transpose() # normalize each GLMC if normed: - out = out.astype(np.float64) + hist = hist.astype(np.float64) for d in range(len(distances)): for a in range(len(angles)): - if np.any(out[:, :, d, a]): - out[:, :, d, a] /= out[:, :, d, a].sum() + if np.any(hist[:, :, d, a]): + hist[:, :, d, a] /= hist[:, :, d, a].sum() - return out + return hist def compute_glcm_prop(P, prop='contrast'): From ab780b1ecd9a59a168aab320477e297a50acf77c Mon Sep 17 00:00:00 2001 From: Neil Yager Date: Mon, 14 Nov 2011 14:32:54 +0000 Subject: [PATCH 07/13] 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 From 0597cec0061a3bba17d86778cfdd70910db1d4f7 Mon Sep 17 00:00:00 2001 From: Neil Yager Date: Mon, 21 Nov 2011 10:06:59 +0000 Subject: [PATCH 08/13] BUG: Fix typo and docstrng --- doc/examples/plot_glcm.py | 2 +- skimage/feature/greycomatrix.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/examples/plot_glcm.py b/doc/examples/plot_glcm.py index e30f10d7..8871a14a 100644 --- a/doc/examples/plot_glcm.py +++ b/doc/examples/plot_glcm.py @@ -90,5 +90,5 @@ plt.ylabel('GLVM Correlation') plt.legend() # display the patches and plot -plt.suptitle('Grey level co-occurance matrix features', fontsize=14) +plt.suptitle('Grey level co-occurrence matrix features', fontsize=14) plt.show() diff --git a/skimage/feature/greycomatrix.py b/skimage/feature/greycomatrix.py index 00929ee8..2b98386d 100644 --- a/skimage/feature/greycomatrix.py +++ b/skimage/feature/greycomatrix.py @@ -156,7 +156,7 @@ def greycoprops(P, prop='contrast'): ... [2, 2, 3, 3]], dtype=np.uint8) >>> g = greycomatrix(image, [1, 2], [0, np.pi/2], levels=4, ... normed=True, symmetric=True) - >>> contrast = compute_glcm_prop(g, 'contrast') + >>> contrast = greycoprops(g, 'contrast') >>> contrast array([[ 0.58333333, 1. ], [ 1.25 , 2.75 ]]) From 830c6b650b8aad4c1a1a79842233515d053ebef3 Mon Sep 17 00:00:00 2001 From: Neil Yager Date: Tue, 22 Nov 2011 11:19:03 +0000 Subject: [PATCH 09/13] ENH: Address division by zero errors --- skimage/feature/greycomatrix.py | 53 +++++++++++++++++------------- skimage/feature/tests/test_glcm.py | 2 +- 2 files changed, 32 insertions(+), 23 deletions(-) diff --git a/skimage/feature/greycomatrix.py b/skimage/feature/greycomatrix.py index 2b98386d..06287fca 100644 --- a/skimage/feature/greycomatrix.py +++ b/skimage/feature/greycomatrix.py @@ -1,5 +1,5 @@ """ -Compute grey level co-occurrence matrices (GLCMs) and associated +Compute grey level co-occurrence matrices (GLCMs) and associated properties to characterize image textures. """ @@ -9,7 +9,7 @@ import skimage.util from _greycomatrix import _glcm_loop -def greycomatrix(image, distances, angles, levels=256, symmetric=False, +def greycomatrix(image, distances, angles, levels=256, symmetric=False, normed=False): """Calculate the grey-level co-occurrence matrix. @@ -102,8 +102,15 @@ def greycomatrix(image, distances, angles, levels=256, symmetric=False, # normalize each GLMC if normed: P = P.astype(np.float64) - P /= np.apply_over_axes(np.sum, P, axes=(0, 1)) - P = np.nan_to_num(P) + glcm_sums = np.apply_over_axes(np.sum, P, axes=(0, 1)) + if np.any(glcm_sums == 0): + # GLCMs are sometimes all zero, so temporarily suppress warning + old_settings = np.seterr(invalid='ignore') + P /= glcm_sums + np.seterr(invalid=old_settings['divide']) + P = np.nan_to_num(P) + else: + P /= glcm_sums return P @@ -190,24 +197,26 @@ def greycoprops(P, prop='contrast'): 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() - 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()) - cov = (g * (diff_i * diff_j)).sum() - if std_i < 1e-15 or std_j < 1e-15: - corr = 1. - else: - corr = cov / (std_i * std_j) - - results[d, a] = corr - - results[d, a] = corr + 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] diff --git a/skimage/feature/tests/test_glcm.py b/skimage/feature/tests/test_glcm.py index 26c4a1f0..3f321e55 100644 --- a/skimage/feature/tests/test_glcm.py +++ b/skimage/feature/tests/test_glcm.py @@ -134,7 +134,7 @@ class TestGLCM(): def test_uniform_properties(self): im = np.ones((4, 4), dtype=np.uint8) - result = greycomatrix(im, [1, 2], [0, np.pi / 2], 4, normed=True, + result = greycomatrix(im, [1, 2, 8], [0, np.pi / 2], 4, normed=True, symmetric=True) for prop in ['contrast', 'dissimilarity', 'homogeneity', 'energy', 'correlation', 'ASM']: From 67b408a3fc0de4365884b10bf0ed20f9efc8f209 Mon Sep 17 00:00:00 2001 From: Neil Yager Date: Tue, 22 Nov 2011 11:55:40 +0000 Subject: [PATCH 10/13] BUG: Fixed warning suppression error --- skimage/feature/greycomatrix.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/feature/greycomatrix.py b/skimage/feature/greycomatrix.py index 06287fca..41bcc3a8 100644 --- a/skimage/feature/greycomatrix.py +++ b/skimage/feature/greycomatrix.py @@ -107,7 +107,7 @@ def greycomatrix(image, distances, angles, levels=256, symmetric=False, # GLCMs are sometimes all zero, so temporarily suppress warning old_settings = np.seterr(invalid='ignore') P /= glcm_sums - np.seterr(invalid=old_settings['divide']) + np.seterr(invalid=old_settings['invalid']) P = np.nan_to_num(P) else: P /= glcm_sums From 7a191d8d20eb2755e7f568988275983ae507ae15 Mon Sep 17 00:00:00 2001 From: Neil Yager Date: Wed, 23 Nov 2011 11:23:03 +0000 Subject: [PATCH 11/13] Long lines, remove doctest, division by zero, data type ranges --- skimage/feature/greycomatrix.py | 37 +++++++++++++++------------------ 1 file changed, 17 insertions(+), 20 deletions(-) diff --git a/skimage/feature/greycomatrix.py b/skimage/feature/greycomatrix.py index 41bcc3a8..907661a2 100644 --- a/skimage/feature/greycomatrix.py +++ b/skimage/feature/greycomatrix.py @@ -18,9 +18,9 @@ def greycomatrix(image, distances, angles, levels=256, symmetric=False, Parameters ---------- - image : ndarray - Input image. The image is converted to the uint8 data type, so - its range of the image is [0, 255]. + image : array_like + 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 @@ -28,7 +28,8 @@ def greycomatrix(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 maximum value is + 256, and the default 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 @@ -80,10 +81,13 @@ def greycomatrix(image, distances, angles, levels=256, symmetric=False, [0, 0, 0, 0]], dtype=uint32) """ - image = np.ascontiguousarray(skimage.util.img_as_ubyte(image)) + + 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 @@ -103,14 +107,8 @@ def greycomatrix(image, distances, angles, levels=256, symmetric=False, if normed: P = P.astype(np.float64) glcm_sums = np.apply_over_axes(np.sum, P, axes=(0, 1)) - if np.any(glcm_sums == 0): - # GLCMs are sometimes all zero, so temporarily suppress warning - old_settings = np.seterr(invalid='ignore') - P /= glcm_sums - np.seterr(invalid=old_settings['invalid']) - P = np.nan_to_num(P) - else: - P /= glcm_sums + glcm_sums[glcm_sums == 0] = 1 + P /= glcm_sums return P @@ -123,11 +121,13 @@ def greycoprops(P, prop='contrast'): 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}\\left|i-j\\right|` + - '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]` + - '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 @@ -138,7 +138,8 @@ def greycoprops(P, prop='contrast'): `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 + prop : {'contrast', 'dissimilarity', 'homogeneity', 'energy', \ + 'correlation', 'ASM'}, optional The property of the GLCM to compute. The default is 'contrast'. Returns @@ -222,7 +223,3 @@ def greycoprops(P, prop='contrast'): results = np.apply_over_axes(np.sum, (P * weights), axes=(0, 1))[0, 0] return results - -if __name__ == "__main__": - import doctest - doctest.testmod() From 7df848623122de6540534f8afaababccbc93b2ac Mon Sep 17 00:00:00 2001 From: Neil Yager Date: Mon, 28 Nov 2011 11:22:31 +0000 Subject: [PATCH 12/13] Update to docstring --- skimage/feature/greycomatrix.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/feature/greycomatrix.py b/skimage/feature/greycomatrix.py index 907661a2..622787f1 100644 --- a/skimage/feature/greycomatrix.py +++ b/skimage/feature/greycomatrix.py @@ -18,7 +18,7 @@ def greycomatrix(image, distances, angles, levels=256, symmetric=False, Parameters ---------- - image : array_like + 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 @@ -29,7 +29,7 @@ def greycomatrix(image, distances, angles, levels=256, symmetric=False, 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, and the default is 256. + 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 From 917d1a80e895ac49cd55d5a27d8ba068a7eb94d8 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Wed, 30 Nov 2011 09:48:01 -0800 Subject: [PATCH 13/13] DOCS: Tighten borders around image in GLCM example. --- doc/examples/plot_glcm.py | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/examples/plot_glcm.py b/doc/examples/plot_glcm.py index 8871a14a..2c848523 100644 --- a/doc/examples/plot_glcm.py +++ b/doc/examples/plot_glcm.py @@ -78,6 +78,7 @@ for (y, x) in sky_locations: plt.xlabel('Original Image') plt.xticks([]) plt.yticks([]) +plt.axis('image') # for each patch, plot (dissimilarity, correlation) plt.subplot(3, 2, 2)