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