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/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..2c848523 --- /dev/null +++ b/doc/examples/plot_glcm.py @@ -0,0 +1,95 @@ +""" +===================== +GLCM Texture Features +===================== + +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 +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. + +""" + +from skimage.feature import greycomatrix, greycoprops +from skimage import data +import matplotlib.pyplot as plt + +PATCH_SIZE = 21 + +# open the camera image +image = data.camera() + +# select some patches from grassy areas of the image +grass_locations = [(474, 291), (440, 433), (466, 18), (462, 236)] +grass_patches = [] +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 +sky_locations = [(54, 48), (21, 233), (90, 380), (195, 330)] +sky_patches = [] +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 = 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 +for i, patch in enumerate(grass_patches): + 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)) + +for i, patch in enumerate(sky_patches): + 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([]) +plt.axis('image') + +# for each patch, plot (dissimilarity, correlation) +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', + label='Sky') +plt.xlabel('GLCM Dissimilarity') +plt.ylabel('GLVM Correlation') +plt.legend() + +# display the patches and plot +plt.suptitle('Grey level co-occurrence matrix features', fontsize=14) +plt.show() diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index eb6faa62..6b3b7014 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 greycomatrix, greycoprops 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 new file mode 100644 index 00000000..622787f1 --- /dev/null +++ b/skimage/feature/greycomatrix.py @@ -0,0 +1,225 @@ +""" +Compute grey level co-occurrence matrices (GLCMs) and associated +properties to characterize image textures. +""" + +import numpy as np +import skimage.util + +from _greycomatrix import _glcm_loop + + +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 + greyscale values at a given offset over an image. + + Parameters + ---------- + 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 + List of pixel pair distance offsets. + angles : array_like + 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 + (typically 256 for an 8-bit image). The maximum value 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 + (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 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 + ------- + 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 + grey-level `i`. If `normed` is `False`, the output is of + type uint32, otherwise it is float64. + + References + ---------- + .. [1] The GLCM Tutorial Home Page, + 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 = greycomatrix(image, [1], [0, np.pi/2], levels=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) + + """ + + 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 + assert angles.ndim == 1 + + P = np.zeros((levels, levels, len(distances), len(angles)), + dtype=np.uint32, order='C') + + # count co-occurences + _glcm_loop(image, distances, angles, levels, P) + + # make each GLMC symmetric + if symmetric: + P += np.transpose(P, (1, 0, 2, 3)) + + # normalize each GLMC + if normed: + P = P.astype(np.float64) + glcm_sums = np.apply_over_axes(np.sum, P, axes=(0, 1)) + glcm_sums[glcm_sums == 0] = 1 + P /= glcm_sums + + return P + + +def greycoprops(P, prop='contrast'): + """Calculate texture properties of a GLCM. + + 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: + + - '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}|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] + + + Parameters + ---------- + 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 + grey-level i. + prop : {'contrast', 'dissimilarity', 'homogeneity', 'energy', \ + 'correlation', 'ASM'}, optional + The property of the GLCM to compute. The default is 'contrast'. + + Returns + ------- + results : 2-D 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 + -------- + 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 = greycomatrix(image, [1, 2], [0, np.pi/2], levels=4, + ... normed=True, symmetric=True) + >>> contrast = greycoprops(g, 'contrast') + >>> contrast + array([[ 0.58333333, 1. ], + [ 1.25 , 2.75 ]]) + + """ + + 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 + + # create weights for specified property + I, J = np.ogrid[0:num_level, 0:num_level] + 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 + 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) + 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] + + return results 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()) + ) diff --git a/skimage/feature/tests/test_glcm.py b/skimage/feature/tests/test_glcm.py new file mode 100644 index 00000000..3f321e55 --- /dev/null +++ b/skimage/feature/tests/test_glcm.py @@ -0,0 +1,144 @@ +import numpy as np +from skimage.feature import greycomatrix, greycoprops + + +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): + 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], + [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): + 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], + [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_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 = greycomatrix(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) + + def test_output_combo(self): + im = np.array([[0], + [1], + [2], + [3]], dtype=np.uint8) + 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) + 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_output_empty(self): + result = greycomatrix(self.image, [10], [0], 4) + np.testing.assert_array_equal(result[:, :, 0, 0], + np.zeros((4, 4), dtype=np.uint32)) + 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_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 = greycomatrix(self.image, [1, 2], [0], 4, + normed=True, symmetric=True) + result = np.round(result, 3) + contrast = greycoprops(result, 'contrast') + np.testing.assert_almost_equal(contrast[0, 0], 0.586) + + def test_dissimilarity(self): + result = greycomatrix(self.image, [1], [0, np.pi / 2], 4, + normed=True, symmetric=True) + result = np.round(result, 3) + dissimilarity = greycoprops(result, 'dissimilarity') + np.testing.assert_almost_equal(dissimilarity[0, 0], 0.418) + + def test_dissimilarity_2(self): + result = greycomatrix(self.image, [1, 3], [np.pi/2], 4, + normed=True, symmetric=True) + result = np.round(result, 3) + dissimilarity = greycoprops(result, 'dissimilarity')[0, 0] + np.testing.assert_almost_equal(dissimilarity, 0.664) + + def test_invalid_property(self): + result = greycomatrix(self.image, [1], [0], 4) + np.testing.assert_raises(ValueError, greycoprops, + result, 'ABC') + + def test_homogeneity(self): + result = greycomatrix(self.image, [1], [0, 6], 4, normed=True, + symmetric=True) + homogeneity = greycoprops(result, 'homogeneity')[0, 0] + np.testing.assert_almost_equal(homogeneity, 0.80833333) + + def test_energy(self): + result = greycomatrix(self.image, [1], [0, 4], 4, normed=True, + symmetric=True) + energy = greycoprops(result, 'energy')[0, 0] + np.testing.assert_almost_equal(energy, 0.38188131) + + def test_correlation(self): + result = greycomatrix(self.image, [1, 2], [0], 4, normed=True, + symmetric=True) + 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 = greycomatrix(im, [1, 2, 8], [0, np.pi / 2], 4, normed=True, + symmetric=True) + for prop in ['contrast', 'dissimilarity', 'homogeneity', + 'energy', 'correlation', 'ASM']: + greycoprops(result, prop) + +if __name__ == '__main__': + np.testing.run_module_suite()