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